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Simulating  saturated/unsaturated  flow  between  wetlands  and  the  contiguous  subsurface 

is  complex,  because  it  involves  the  modeling  of  wetland  storage  changes  that  in  turn  cause 

changes  in  pertinent  hydrologic  boundary  conditions.  Existing  groundwater  models  are 

incapable  of  simulating  the  flow  between  a  wetland  and  an  aquifer,  if  small  changes  in 

wetland  storage  cause  significant  vertical  and  lateral  movements  of  the  wetland  boundary. 

In  addition  to  the  simulation  complexities  introduced  by  a  change  in  wetland  geometry  and 

storage,  periodic  and  nonperiodic  fluctuations  of  the  wetland  stage  induce  corresponding 

elevation  changes  in  the  phreatic  surface  of  the  contiguous  aquifer.  As  these  induced  waves 

propagate  through  an  aquifer,  friction  causes  a  loss  of  energy,  which  is  manifested  as 

spatially  dampened  potentiometric  head  perturbations  along  the  inland  direction.  An 

analytical  model  was  developed  that  describes  subsurface  flows  around  a  wetland  induced 

by  any  nonperiodic  fluctuations  in  surface-water  stage,  such  as  a  flood  wave.  Not  considered, 
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however,  are  any  changes  in  wetland  storage,  or  wetland  geometry  as  a  function  of  surface 
water  stage.  The  analytical  model  was  validated  using  a  finite-difference  numerical 
groundwater  model  that  is  based  on  the  governing  equation  expressed  in  radial  coordinates. 
Comparisons  of  analytical  and  numerical  results  show  excellent  agreement. 

Two  numerical  models  were  specifically  designed  to  simulate  wetland-aquifer 
interactions.  The  first  model  is  a  saturated  groundwater  flow  model  that  incorporates 
adaptive  simulation  technologies  to  permit  real-time  simulation  of  moving  wetland 
boundaries  and  their  associated  local  influence  on  the  phreatic  surface.  This  model  is 
numerically  efficient  and  uses  deformable  hexahedral  finite  elements  to  trace  phreatic  surface 
changes  in  real-time.  The  second  model  also  uses  deformable  hexahedral  finite  elements; 
however,  this  model  simulates  variably  saturated  groundwater  flow.  For  both  models,  the 
numerical  elements  deform  as  required  to  characterize  the  horizontal  and  vertical  extent  of 
the  moving  wetland  boundary  (which  serves  to  define  a  location  within  the  porous  system, 
where  the  pressure  is  known).  Simulation  results  from  both  models  were  applied  to  a  field 
site  where  water  was  rapidly  withdrawn  from  an  isolated  wetland  to  observe  system  response 
and  evaluate  the  wetland  aquifer  hydraulic  connection. 
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CHAPTER  1 
INTRODUCTION 

1.1  Wetland 

Wetlands  are  transitional  lands  between  terrestrial  and  aquatic  systems  where  the 
water  table  is  usually  at  or  near  the  surface,  or  the  land  is  covered  by  shallow  water 
(Cowardin  et  al.  1979).  Wetlands  are  also  called  marsh,  fen,  peatland  or  water  and  can  be 
natural  or  artificial  and  permanent  or  temporary.  The  water  can  be  static  or  flowing  and  fresh 
or  brackish,  or  it  can  include  areas  of  marine  water  the  depth  of  which  at  low  tide  does  not 
exceed  6  meters  (Ramsar  Convention  1971,  Middleton  1998). 

Wetlands  perform  various  environmental  functions  that  include  increasing  flood 
mitigation,  enhancing  storm  abatement  (i.e.,  in  the  case  of  coastal  wetlands),  adding 
groundwater  recharge,  providing  water  quality  treatment,  serving  as  wildlife  habitats,  and 
performing  aesthetic  functions  (Mitsch  and  Gosselink  1993).  These  functional  attributes  can 
be  classified  into  three  basic  types:  hydrologic,  water  quality,  and  life  support.  Measures  of 
the  life-support  function  must  integrate  several  wetland  attributes,  including  water  quality, 
hydroperiod  and  hydrodynamics,  habitat  structure,  food  chain  support,  biogeographic  setting, 
nutrient  status,  corridors  for  migration,  and  primary  production  (Preston  and  Bedford  1988). 

In  recent  years,  increasing  interest  has  focused  on  the  study  and  protection  of 
wetlands.  Development,  drainage,  population  pressure,  and  pollution  have  destroyed 
extensive  areas  of  these  sensitive  ecosystems.  Another  increasingly  significance  to  wetlands 
is  municipal  well  fields  supplying  water  to  large  urban  areas.  These  well  fields  affect 
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2 
wetlands  by  reducing  surficial  aquifer  levels  over  large  areas  that  encompass  isolated  and 
interconnected  wetlands  (Sonenshein  and  Hofstetter  1990).  As  a  result  of  depressing  the 
water  table,  historically  wet  ecosystems  have  been  drained  (Switt  et  al.  1998). 

1.2  Wetland  Water  Budget 
Hydrologic  conditions  are  extremely  important  for  the  maintenance  of  a  wetland' s 
structure  and  function.  A  hydrologic  budget  is  often  used  to  quantify  the  volume  and  rate  of 
water  flowing  into  and  out  of  a  wetland.  The  primary  hydrologic  inflows  include 
precipitation,  streams,  groundwater,  and,  in  coastal  wetlands,  tides  (Mitsch  and  Gosselink 
1986,  Rosenberry  2000).  The  wetland  water  budget  equation  as  defined  by  Mitsch  and 
Gosselink  (1993)  is: 


AF/Af  =  Pn+St  +  GrET-So  -  G±T 


(1-2) 


where 

AV  =  change  in  the  volume  of  water  storage  in  wetland  [L3] 

At  =  change  in  time  [T] 

Pn  =  net  precipitation  ( total  precipitation  minus  intercepted  precipitation  (I) )  [  L3/T] 

Sj  =  surface  inflows,  including  flooding  stream  [  L3/T] 

Gj  =  groundwater  inflow  [L3/T] 

ET  =  evapotranspiration  [L3/T] 

S  0  =  surface  outflows  [L3/T] 

G0  =  groundwater  outflows  [L3/T] 

T  =  tidal  inflow  (+)  or  outflow  (-)  [L3/T] 
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Terms  in  the  wetland  water  budget  vary  in  importance  depending  on  the  wetland  type, 
and  not  all  terms  are  required  to  characterize  a  specific  type  of  wetland.  In  most  wetland 
systems,  precipitation  and  evapotranspiration  are  major  components  of  the  water  budget 
(Winter  and  Woo  1990).  Groundwater  can  also  have  a  significant  influence  on  some  wetland 
systems,  but  not  always. 

Three  of  the  most  common  relationships  that  characterize  a  wetland  are  groundwater 
discharge,  groundwater  inflow  and  outflow,  and  groundwater  recharge.  Groundwater  inflow 
occurs  when  the  contiguous  wetland  surface  water  (or  groundwater)  stage  is  lower  than  the 
phreatic  surface  of  the  surrounding  aquifer.  When  the  surface  water  stage  or  groundwater 
levels  inside  a  wetland  are  higher  than  the  phreatic  surface  in  the  surrounding  contiguous 
aquifer,  groundwater  and  surface  water  will  flow  out  of  wetland.  When  a  wetland  is  situated 
well  above  the  groundwater,  the  wetland  is  identified  as  being  perched.  This  type  of 
hydrologic  system,  referred  to  as  a  surface  water  depression  wetland  by  Novitzki  (1979), 
loses  water  only  through  infiltration  into  the  ground  and  through  evapotranspiration. 

The  rate  of  water  exchange  between  a  wetland  and  local  aquifer  is  a  function  of  the 
hydraulic  gradient  between  the  two  hydrologic  subsystems.  The  discharge,  in  accordance 
with  Darcy's  law  (Truesdell  1995),  is  proportional  to  the  hydraulic  gradient,  the  hydraulic 
conductivity,  and  the  cross-sectional  area  through  which  the  flow  occurs.  Periodic  or 
nonperiodic  fluctuations  of  the  wetland  stage  in  a  flow-through  wetland  will  induce 
fluctuations  in  the  elevation  of  the  phreatic  surface  in  the  contiguous  aquifer.  As  these 
induced  waves  propagate  through  an  aquifer,  friction  causes  a  loss  of  energy,  which  is 
manifested  as  spatially  dampened  head  perturbations.  This  damping  effect  progressively 
diminishes  the  amplitude  of  piezometric  head  along  the  inland  direction  (Carr  and  Van  Der 
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Kamp  1969).  Many  studies  have  examined  the  problem  through  one-dimensional  analytical 
solutions  (Cooper  1963,  Gregg  1966,  Nielsen  1990,  Serfes  1992,  Li,  Barry  et  al.  2000, 2001, 
Boufadel  and  Peridier  2002,  Seerano  2003,  Swamee  and  Singh  2003).  A  two-dimensional 
analytical  solution  was  obtained  by  considering  a  shift  in  phase  and  a  change  in  amplitude 
along  a  straight  coastline  (Sun  1997).  Three-dimensional  flow  regimes  near  a  shallow 
circular  lake  was  studied  using  numerical  and  analytical  mathematical  models  of  aquifer  flow 
(Townley  and  Trefry  2000). 

To  further  the  understanding  of  the  interactions  between  lakes  or  wetlands  and 
groundwater,  without  the  restrictive  assumptions  typically  associated  with  analytical 
solutions,  numerical  models  are  usually  used.  For  example,  it  is  possible  to  obtain  numerical 
solutions  for  the  case  of  anisotropic  and  nonhomogeneous  conditions.  The  accuracy  of 
solutions  obtained  by  numerical  methods  can  be  excellent;  however,  it  depends  on  several 
factors,  including  the  type  of  numerical  method  used,  the  complexity  of  boundary  and  initial 
conditions,  and  the  computational  precision  of  the  computer  used  to  implement  the  method. 
Regardless  of  these  advantages,  few  investigators  have  chosen  to  focus  on  the  numerical 
simulation  of  interactions  between  wetlands  and  groundwater  (  McDonald  and  Harbaugh 
1988,  Schaffranek  1996,  Lee  2000). 

1.3  Goals  and  Objectives 

Analytical  models  to  simulate  the  complex  groundwater  surface  water  interactions 
that  occur  with  isolated  and  flow-through  wetlands  would  be  valuable  tools  for  developing 
hydrologic  budgets  for  wetlands.  Given  that  most  surface  water  fluctuations  are  nonperiodic, 
there  is  a  need  to  develop  analytical  groundwater  models  for  simulating  groundwater  table 
fluctuations  induced  by  both  periodic  and  nonperiodic  surface  water  flood  waves.  However, 
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with  the  exception  of  the  one-dimensional  study  by  Cooper  and  Rorabaugh  (1963),  most 

investigations  addressed  surface-water  fluctuations  that  were  periodic.  Thus,  the  first  goal 
of  this  dissertation  was  to  develop  an  analytical  method  for  simulating  groundwater  table 
fluctuations  induced  by  both  periodic  and  nonperiodic  surface-water  flood  waves.  This  goal 
was  achieved  by  meeting  the  single  objective  of  developing  an  analytical  model  using 
Fourier  transforms. 

A  limited  number  of  three-dimensional  numerical  groundwater  models  have  been 
developed  to  simulate  groundwater  wetland  interactions;  however,  these  models  are  generally 
inadequate  tools  for  simulating  the  complex  interactions  between  a  wetland  and  a  contiguous 
unconfined  aquifer.  For  example,  Walser,  Annable  and  Wise  (1998)  used  FEMWATER 
(Yeh,  1992),  a  variably  saturated  flow  model,  to  simulate  the  exchange  of  water  between  a 
wetland  and  an  aquifer.  McDonald  and  Harbaugh  (1988)  used  MODFLOW  and  developed 
a  three-dimensional  finite  difference  groundwater  model  and  studied  the  dynamic  coupling 
between  wetlands  and  the  subsurface.  Lee  (2000)  used  the  HYDRUS-2D  flow  model  and 
simulated  vertical  two-dimensional  groundwater  flow  through  Lake  Barco  in  Florida.  Neither 
FEMWATER  (1996)  nor  MODFLOW  can  simulate  or  incorporate  a  moving  real-time 
wetland  boundary  caused  by  a  change  in  surface  water  stage.  Furthermore,  as  a  saturated 
three-dimensional  finite-difference  model,  MODFLOW  can  not  simulate  fluctuations  of  the 
phreatic  surface  as  a  moving  transient  boundary  as  required  to  define  fully  the  three- 
dimensional  flow.  Finally,  the  water  stage  of  a  wetland  itself  is  a  function  of  the  transient 
exchange  of  water  between  a  surface  and  subsurface  systems;  which  is  not  fully  modeled  in 
any  existing  groundwater  flow  codes.  Thus,  to  simulate  groundwater  surface  water 
interaction  without  the  above  identified  shortcomings  of  existing  models,  a  new  model  is 


needed. 

The  second  goal  of  the  dissertation  was  to  develop  two  numerical  models  capable 
of  simulating  1)  a  transient  and  moving  wetland  boundary,  2)  three-dimensional  flow,  3)  the 
movement  of  the  phreatic  surface,  and  4)  the  coupling  that  occurs  between  the  transient 
surface-water  stage  and  the  phreatic  surface.  To  achieve  this  goal,  two  objectives  were 
pursued.  The  first  objective  was  to  develop  two  numerical  groundwater  models  using  finite 
elements.  The  second  objective  was  to  apply  both  numerical  models  to  simulate  groundwater 
and  surface  water  fluctuations  over  an  isolated  wetland  in  southern  Florida.  Of  the  two 
models  developed,  the  first  is  a  three-dimensional  saturated  groundwater  flow  model  that 
employs  deformable  hexahedral  elements.  The  finite-element  mesh  deforms  in  space  and 
time  to  trace  the  location  of  the  phreatic  surface  of  an  unconfined  aquifer  and  the  stage  of  the 
contiguous  surface  water  body.  Based  on  the  above  model  and  the  adjoint  state  method,  an 
inverse  model  was  developed  to  estimate  the  hydraulic  conductivity  of  peat  layers  which 
typically  exist  as  a  wetland  bottom  that  serves  to  isolate  surface  wetlands  from  subsurface 
groundwater  systems.  The  second  numerical  model,  a  three-dimensional  variable  saturated 
flow  model,  was  also  developed  using  deformable  mesh.  This  model  was  used  to  investigate 
prediction  differences  between  the  saturated  flow  code  and  the  variable  saturated  model. 
From  the  groundwater  flow  pattern  given  by  both  models,  an  improved  understanding  was 
obtained  for  the  interactions  occurring  between  surface  waters  and  groundwaters.  In  addition, 
an  improved  understanding  was  obtained  of  how  short-term  pumping  from  a  phreatic  aquifer 
affects  wetlands  and  lake  water  stages. 


CHAPTER  2 
ANALYTICAL  GROUNDWATER  MODEL  RESPONSE  TO  THE  PERIODIC 
OR  NON-PERIODIC  WATER  STAGE  FLUCTUATIONS  IN  A  CIRCULAR  WETLAND 


2.1  Introduction 

Periodic  or  nonperiodic  fluctuations  of  the  surface  water  stage  in  wetlands  or  lakes 
will  induce  fluctuations  in  the  elevation  of  the  phreatic  surface  of  contiguous  unconfined 
aquifers.  As  these  induced  waves  propagate  through  the  aquifer,  friction  causes  a  loss  of 
energy,  which  is  manifested  as  spatially  dampened  head  perturbations.  This  damping  effect 
progressively  diminishes  the  amplitude  of  fluctuations  in  the  piezometric  head  along  the 
inland  direction  (Carr  and  Van  Kamp  1969).  Many  studies  have  examined  the  problem 
through  one-dimensional  analytical  solutions  (e.g.,Cooper  and  Rorabaugh  1963,  Gregg  1966, 
Nielsen  1990,  and  Serfes  1992,  Li,  Barry  et  al.  2000,  2001,  Boufadel  and  Peridier  2002, 
Seerano  2003,  Swamee  and  Singh  2003).  Cooper  and  Rorabaugh  (1963)  obtained  the 
analytical  solution  for  groundwater  movement  and  bank  storage  due  to  flood  stages  in  surface 
streams.  Gregg  (1966)  found  a  one-dimensional  relation  between  tidal  fluctuations  and 
groundwater  table  perturbations.  Nielsen  (1990),  Serfes  (1992)  and  Carr  (1969)  also 
published  papers  investigating  the  similar  problems.  A  two-dimensional  analytical  solution 
was  obtained  by  considering  a  phase  shift  in  water-table  perturbations  and  the  associated 
change  of  amplitude  along  a  straight  coastline  (Sun  1997).  An  analytical  solution 
representing  vertical  two-dimensional  groundwater-surface  water  interaction  was  obtained 
by  Anderson  (2003).    Three-dimensional  flow  regimes  near  a  shallow  circular  lake  was 
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studied  using  numerical  and  analytical  mathematical  models  of  aquifer  flow  (Townley  and 
Trefry  2000,  Smith  and  Townley  2002). 

With  the  exception  of  Cooper  and  Rorabaugh  (1963),  all  of  the  above  studies 
addressed  periodic  surface-water  fluctuations,  in  spite  of  the  fact  that  most  stage  variations 
are  not  periodic  i.e.,  consider,  for  example,  the  pervasive  nonperiodic  flooding  of  rivers, 
lakes,  and  wetlands. 

Cooper's  (1963)  nonperiodic  analytical  model  for  a  river/groundwater  system  was 
a  one-dimensional  solution  along  an  axis  perpendicular  to  the  river.  However,  this  analytical 
model  is  not  suitable  for  simulating  the  nonperiodic  interactions  between  a  wetland  and  an 
aquifer.  In  this  chapter,  two  analytical  groundwater  models  are  developed  to  describe 
transient  groundwater  levels  affected  by  periodic  and  nonperiodic  fluctuations  of  the  surface- 
water  stage  in  a  wetland  or  a  lake. 

2.2  Problem  Formulation:  Groundwater  Flow  Around  a  Circular  Wetland 

Considering  a  circular  wetland  that  is  hydraulically  contiguous  to  a  surfacial  aquifer 
and  ignoring  recharge  from  rainfall,  the  governing  equation  for  two-dimensional  groundwater 
flow  in  a  homogeneous  isotropic  aquifer  can  be  expressed  in  radial  coordinates  as: 

32h     1  dh    Sdh 

dr2    r  dr    T  at 

where  h  is  the  piezometric  head,  r  is  the  radial  coordinate,  t  is  time,  S  is  the  specific  storage 
and  Tis  the  transmissivity  in  the  aquifer.  When  a  thick  unconfined  aquifer  is  considered,  the 
gradient  and  small  fluctuations  relative  to  saturated  thickness  are  assumed  to  be  small,  Tis 
approximately  equal  to  Khm  (Figure  2-1)  in  which  K  is  the  hydraulic  conductivity  and  hno  is 


the  average  unconfined  aquifer  thickness  equal  to  the  regional  average  water-table  elevation 
above  an  underlying  impermeable  confining  unit. 

The  boundary  condition  along  the  shore  line  of  the  wetland  is  Eq.  2-2,  where  r  0  is 
the  radius  of  the  circular  wetland.  It  is  also  assumed  in  Eq.  2-3  that  as  r-°°,  the  groundwater 
table  elevation  approaches  the  regional  average  water  table  elevation,  h^. 


B.C:  h(r,t)\  r=ro  =  h„  +f(r0  ,t)  @  r=  r0 


(2-2) 


h(r,t)\r=„  =h(, 


@  r= 


(2-3) 


wetland  or  lake 


//////////////////////////////// 


Figure  2-1.  Setup  of  the  groundwater  problem  around  a  circular  wetland 
2.3  Analytical  Solution  for  the  Periodic  Boundary  Condition 

Since  the  boundary  condition  is  a  periodic  function  of  time,  substituting  for  f(rott), 
Equation  2-4  is  obtained: 

B.C:  h(r,t)\  r=ro  =  hm  +  h0eia'  @  r=  r0  (2-4) 

h(r,t)\r=„=h00  @r=°o  (2-5) 

where  h0  is  the  half  height  of  the  periodic  constituent  and  w  is  the  angular  frequency.  No 

initial  condition  is  necessary,  since  here  the  steady-state  oscillatory  solution  is  sought  as  t-°°. 

Because  of  the  oscillatory  boundary  condition,  a  solution  is  sought  in  the  form 
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h(r,t)=hon+hog(r)e^ 
Substitution  of  Eq.  2-6  into  Eq.  2-1,  Eq.  2-2,  and  Eq.  2-3  yields 


v= 


(2-6) 


&g+Ldg_Sa       Q  (2.7) 

dr2    r  dr      T 


where  the  boundary  conditions  are: 

g(r)=l        @  r=r„  (2-8) 

g(r)=0        @  r=°°  (2-9) 

Equation  2-7  can  now  be  rewritten  to  obtain: 


dv2      dv       *  (2-10) 


N^('"1)T  (2-1D 


The  solution  to  Equation  2-10  subject  to  the  boundary  conditions  in  Equation  2-8  and  2-9 
depends  on  the  value  of  w.  Considering  the  boundary  condition  along  the  shoreline, 
groundwater  around  a  circular  water  body  will  have  the  following  forms. 
When  a)  is  greater  than  0: 


r/(1Yv) 

g{r)=CQ(Jo+iYo)=C0H*Xv)  =  -jLi!  (2-12) 
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and  when  oo  is  less  than  0: 


(2)  Hf\U) 

g(r)=C0{J0-iYo)=C0H?\u)= 


Hf\u0) 


(2-13) 


where  J0  is  the  Bessel  function  of  the  first  kind  order  0;  Y0  is  the  Bessel  function  of  the 
second  kind  order  0;  H0(t)  is  the  Hankel  function  order  0;  and  C0  is  a  coefficient  (Lebedev 
1965).  Substitution  of  Eq.  2-12  and  Eq.  2-13  into  Eq.  2-6  gives  the  following  solution  for  the 
groundwater  flow  around  a  circular  lake  affected  by  periodic  function: 


(2-14) 


for  G)  >  0  and: 


Hf\u) 


0  0  O  f)\ 


(2-15) 


for  oo<  0. 


The  asymptotic  approximations  of  the  functions  H0(1,(z)  and  H0(2,(z)  are: 
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(2izyk+0(\z\-"-1) 


(2-17) 
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2.4  Analytical  Solution  for  the  Non-Periodic  Boundary  Condition 

A  flood  wave  will  typically  induce  a  nonperiodic  variation  in  the  wetland  stage.  Here, 
f(r0,t),  which  is  a  nonperiodic  function  of  time,  is  assumed  to  describe  transient  variations 
in  the  wetland  stage.  Using  the  Fourier  transformation,  f(r0  ,t)  can  be  represented  as: 


ArJ)=j-\F{<*)emd«> 


(2-18) 


where  F(go)  is  the  Fourier  coefficient,  and  F(oo)doo  is  the  amplitude  of  oscillation  for  the 
component  with  a  frequency  in  the  range  of  (co,oo+d(o).  For  this  F(co),  the  corresponding 
fluctuations  in  the  groundwater  table  are  described  by: 


#f(v) 

"0()(v0) 


(G>>0) 


(2-19) 
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or 


(  oo<  0  )         (2-20) 


as  given  byEq.  2-14  andEq.  2-15.  The  inverse  Fourier  transformation  of  Equation  2-19  and 
Eq.  2-20  leads  to: 


Hr,t)=h00  + 


2% 


F(o>)     °  K     •<?""•<&)  +    F((o)     °  — -e 


""•rfco 


=*"+£ 


IXr0,T)e-"'rfu 


o)/i 


J(0  + 


'^.^-"'rfc 


••« 


Ull 


J(0 


--Ko+^lfcorf 


\    0 


■e"{'-')i-d<x>  + 


Hf\u) 


e  U('-T)/-J(0 


dx 


(2-21) 


From  the  above  analytical  model,  it  is  clear  that  the  unsteady  groundwater  table  around  a 
circular  wetland  depends  not  only  on  the  instantaneous  values  of  the  wetland  water  stage  but 
also  on  the  accumulative  effect  represented  by  the  integration.  However,  physically  it  is  not 
possible  for  the  wetland  water  stage  at  a  later  time  to  affect  the  groundwater  table  at  an 
earlier  time.  Hence,  the  phreatic  surface  depends  on  its  current  stage  and  historical  stage.  The 
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final  analytical  solution  to  the  nonperiodic  boundary  condition  2-2  is  then 


Kr,t)-hoo 


2n 


W) 


Hf\u0) 


e^~z)i-du  + 


e  u('-T)'^co 


dx 


(2-22) 


2.5  Application 
If  it  is  assumed  that  wetland  stage  fluctuations  can  be  described  by  a  nonperiodic 
Gaussian  function  of  time,  then: 


J2-KO 


2o2 


(2-23) 


The  analytical  solution  obtained  when  Eq.  2-23  is  combined  with  Eq.  2-21  is: 


2*  »  v/2™ 


(*-'„)' 


2cr 


Hf\u) 


•ew('-T)'-Jo)  + 


g0(1)(v) 

^U(v0) 


•e  <j('-t)'-Jq 


rft 


(2-24) 
The  transient  changes  of  groundwater  table  described  by  Eq.  2-24  at  four  different 
distances  from  the  wetland  are  shown  in  Figure  2-2.  The  groundwater  distributions  along  the 
radius  at  four  different  times  are  shown  in  Figure  2-3.  From  both  figures,  it  is  evident  that 
groundwater  fluctuations  are  greatest  near  the  wetland.  From  Figure  2-3,  it  can  be  seen  that 
the  variation  of  the  groundwater  table  over  time  is  no  longer  symmetrical  although  the 
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variation  of  the  wetland  water  stage  is  a  Gaussian  distribution  which  is  symmetrical  over 
time. 
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Figure  2-2.  Transient  groundwater  level  changes  at  radial  distances  1,2,3,  and  4  meters 

away  from  the  wetland  shoreline. 
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Figure  2-3.  Groundwater  distribution  over  the  distance  to 
the  wetland/lake  shoreline  at  6,7,8,9,  and  10  h. 
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Figure  2-4  shows  the  groundwater  table  elevations  around  a  circular  wetland  after  3 
hours.  The  radius  of  the  wetland  is  50  meters.  All  other  parameters  are  listed  in  Table  2-1. 
Theoretically,  when  0)  is  very  large,  it  means  that  the  frequency  of  the  corresponding 
boundary  condition  is  very  large  and  the  period  is  very  short.  When  the  absolute  value  of  oa 
is  very  large,  its  effect  on  the  groundwater  around  the  wetland  is  very  small.  This  is  because 
the  period  of  the  wave  component  is  too  small  to  have  any  impact  on  the  groundwater  around 
the  wetland.  That  is  to  say,  when  the  integration  of  co  is  taken,  a  certain  negative  number 
could  be  used  to  replace  the  negative  infinite  number,  and  a  certain  positive  number  could 
also  be  used  to  replace  the  positive  infinite  number.  Accurate  model  results  are  obtained 
when  w  takes  the  range  of  [-5,  5].  Thus,  it  may  be  assumed  that  the  boundary  condition  has 
no  effect  on  the  ground  water  table  if  the  frequency  of  the  boundary  water  level  is  very  large. 


Y  (M) 


Figure  2-4.  Three-dimensional  view  of  the  groundwater  table 
around  a  circular  wetland  at  t=3  h. 
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2.6  Comparison  of  Analytical  Solution  and  Numerical  Solutions 

A  comparison  was  made  between  the  analytical  and  the  numerical  solutions  to  the 
radial  flow  problem  specified  through  Eq.  2-1-Eq.  2-3.  An  implicit  finite-difference  scheme 
was  used  to  represent  the  radial  form  of  the  groundwater  flow  equation  Eq.  2-1  as: 


i  n+l    ~,  n+1     »  n  +  \  ,  n+1     .  n  +  1  i"  +  l     1  " 

+ = 


(Ar)2 


2Ar 


Ar 


(2-25) 


Here  i  is  the  reference  integer  spatial  interval;  n  is  the  current  integer  time  level;  At  is  the 
time  step;  and  Ar  is  the  radial  step. 

A  circular  wetland  which  has  a  radius  of  50  meters  was  assumed.  All  other  system 
parameters  are  given  in  the  Table  2- 1 .  By  considering  the  different  boundary  conditions  (Eq. 
2-3  andEq.  2-4),  the  corresponding  numerical  results  are  found  [Figure  2-5].  Analytical  and 
numerical  results  illustrate  time  variations  in  groundwater  levels.  Figure  2-5  shows  that  the 
analytical  solution  and  numerical  solution  compare  very  well. 

Table  2-1.  Parameters  and  their  values 


Parameter 

Value 

Parameter 

Value 

HO0  (m) 

100.0 

T    (m2/h) 

0.635 

Ar   (m) 

0.25 

S 

0.2 

At    (hr) 

0.02 

t0     (h) 

5.0 

k   (m/hr) 

0.00635 

o 

0.4 
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Figure  2-5.  Comparison  of  analytical  solution  and  numerical  solution 
at  2  and  5  meters  from  the  wetland  shoreline. 


2.7  Summary 

The  analytical  model  presented  in  this  chapter  can  be  used  to  simulate  aquifer  surface 
fluctuations  caused  by  the  transient  wetland  stages  regardless  of  the  temporal  function 
assumed  to  describe  variations  of  water  stage.  In  other  words,  the  model  can  deal  with  any 
general  transient  of  the  boundary  function.  The  analytical  model  was  validated  using  a  finite- 
difference  numerical  groundwater  model  for  a  homogenous  aquifer.  This  first  numerical 
method  is  not  recommended  for  simulations  involving  heterogenous  aquifers.  Rather,  a  new 
method  is  presented  for  simulating  wetland/aquifer  interactions  in  such  aquifers,  which 
requires  a  simple  transformation  of  the  equations.  The  excellent  agreement  between  two 
models  shows  that  the  analytical  solution  presented  is  valid.  This  model  is  also  useful  in  a 
situation  where  groundwater  and  the  surface  water  components  are  being  considered  in  a 
system-wide  water  balance. 


CHAPTER  3  TRANSIENT  THREE-DIMENSIONAL  DEFORMABLE 
FINITE-ELEMENT  SATURATED  GROUNDWATER  MODEL 


3.1  Introduction 

Recent  research  has  focused  on  the  details  of  various  wetland  functions  such  as 
water-quality  treatment,  flood-wave  suppression,  and  groundwater/surface  water  interactions 
(Wise  et  al.  2000,  Krabbenhoft  and  Webster  1995,  Anderson  and  Cheng  1993,  Ruddy  and 
Williams  1991,  and  Mills  and  Zwarich  1986).  For  example,  Ruddy  and  Williams  (1991) 
investigated  the  hydrologic  relation  between  the  South  Fork  Williams  Fork  stream  in  the 
Colorado  River  basin  and  four  adjacent  subalpine  wetlands,  and  they  also  evaluated  the 
potential  effects  on  these  wetlands  when  stream  flow  is  diverted.  They  concluded  that  stream 
flow  was  the  primary  source  of  water  for  two  of  the  four  wetlands,  whereas  precipitation, 
snowmelt,  valley  side-slope  flow  (including  overland  and  small-channel  flow),  and 
groundwater  inflow  served  to  define  the  transient  hydrological  behavior  of  all  wetlands. 

The  above  studies  showed  that  wetland  functions  (i.e.,  water  quality  treatment, 
floodwave  suppression,  etc.)  are  correlated  with  recharge  and  discharge  areas,  ground-water- 
flow  paths  and  rates,  and  soil-profile  development.  The  location  of  recharge  and  discharge 
areas  can  change  in  response  to  seasonal  or  weather-related  factors  and  can  affect  the 
groundwater  flow  direction.  As  a  result,  groundwater  and  surface  water  quality  and  the 
presence  of  various  wetland-plant  communities  may  all  vary  in  response  to  these  changing 
conditions. 
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Analytical  models  such  as  the  one  presented  in  Chapter  2  can  be  used  to  simulate 
deep  wetlands  or  lakes.  However,  the  primary  limitation  of  analytical  methods  is  that 
solutions  can  only  be  obtained  by  imposing  severely  restrictive  assumptions  on  aquifer 
properties  and  boundary  and  initial  conditions.  For  example,  an  assumption  commonly  made 
to  obtain  analytical  solutions  is  that  the  aquifer  is  isotropic  and  homogeneous  for  hydraulic 
conductivity.  These  assumptions,  however,  are  not  valid  and  often  not  adequate  to  describe 
groundwater  flow  or  solute  transport  in  most  field  situations. 

Numerical  methods  have  the  primary  advantage  over  analytical  approaches  in  that 
they  do  not  require  such  restrictive  assumptions.  For  example,  it  is  possible  to  obtain 
numerical  solutions  for  the  case  of  anisotropic  and  nonhomogeneous  conditions.  The 
accuracy  of  solutions  obtained  by  numerical  methods  can  be  excellent.  However,  this 
accuracy  depends  on  several  factors,  including  the  type  of  numerical  method  used,  the 
complexity  of  boundary  and  initial  conditions,  and  the  computational  precision  of  the 
computer  used  to  implement  the  method.  Regardless  of  these  advantage,  only  a  few 
investigators  have  chosen  to  focus  on  the  numerical  simulation  of  interactions  between 
wetlands  and  groundwater  (e.g.,  Schaffranek  1996  and  McDonald  and  Harbaugh  1988,  Lee 
2000,  Rosenberry  2000).  In  this  chapter,  a  general  discussion  of  the 
advantages/disadvantages  of  two  numerical  methods  commonly  used  to  model  hydrologic 
systems  is  briefly  discussed.  Next,  the  mathematical  development  of  a  three-dimensional 
deformable  finite-element  saturated  groundwater  model  is  presented.  This  model  and  the 
associated  code  are  then  validated  by  comparing  model  predictions  with  analytical  modeling 
results.  A  more  complicated  application  of  the  three-dimensional  model  is  presented  to 
illustrate  the  type  of  information  provided  through  the  numerical  approach  taken.  Finally,  the 
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model  is  used  to  simulate  an  actual  wetland  using  field  data  collected  when  the  hydrologic 

system  was  responding  to  induced  stress. 

3.2  Ground-Water  Flow  Models  Based  on  Finite-Difference 
and  Finite-Flement  Approximation 

Several  types  of  numerical  methods  have  been  used  to  solve  groundwater  flow 

problems,  the  two  principal  ones  being  the  finite-difference  method  and  the  finite-element 

method.  The  finite-difference  method  was  initially  applied  to  the  flow  of  fluids  in  petroleum 

reservoirs,  and  it  wasn't  until  the  mid-1960's  that  problems  of  groundwater  flow  and  solute 

transport  were  examined.  The  method  has  a  number  of  advantages  that  contribute  to  its 

continued  widespread  use  and  popularity:  (1)  for  simple  problems  such  as  one-dimensional, 

steady-state  groundwater  flow  in  an  isotropic  and  homogeneous  aquifer,  mathematical 

formulation  and  computer  implementation  are  easily  understood;  and  (2)  the  accuracy  of 

solutions  to  steady-state  and  transient  groundwater  flow  problems  is  generally  quite  good. 

The  finite-difference  method  also  has  disadvantages:  (1)  the  method  works  best  for 

rectangular  or  prismatic  aquifers  of  uniform  composition  (for  example,  it  is  difficult  to 

incorporate  irregular  or  curved  aquifer  boundaries,  anisotropic  and  heterogeneous  aquifer 

properties,  or  sloping  soil  and  rock  layers  into  the  numerical  model  without  introducing 

numerous  mathematical  and  computer  programming  complexities);  and  (2)  the  accuracy  of 

solutions  to  solute  transport  problems  is  less  than  can  be  obtained  through  application  of  the 

finite-element  method.  Table  3.1  lists  selected  early  references  where  finite-difference 

methods  were  applied  to  model  groundwater  flow.  More  recent  applications  of  hydrologic 

modeling  using  finite-difference  for  simulating  of  surface  water  /groundwater  interactions 

include  studies  by  Mc  Donald  and  Harbaugh  (1998),  Winter  (1983, 1986,  and  1989),  Sacks 
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(1992)  and  Krabbenhoft  and  Anderson  (1990). 

Table  3-1.  Selected  early  references  for  the  finite-difference  method. 


Topic 

Reference 

Early  Developments  in  Petroleum 
Reservoir  Modeling 

Bruce  et  al.  (1953),  Peace  and 
Rachford  (1962) 

Saturated  Groundwater  Flow 

Remson  et  al.  (1965),  Freeze  and 
Whitherspoon  (1966),  Pinder  and 
Bredehoeft  (1968). 

Unsaturated  Groundwater  Flow 

Philip  (1957),  Ashroft  et  al.  (1962), 
Freeze  (1971),  Brutsaert  (1973), 
Healy  (1990),  Kirkland  et  al.  (1992). 

Solute  Transport 

Stone  and  Brian  (1963),  Oster  et  al. 
(1970),  Tanji  et  al.  (1967),  Wierenga 
(1977) 

Application  to  Field  Problems 

Orlob  and  Woods  (1967),  Gambolati 
et  al.  (1973),  Fleck  and  McDonald 
(1978),  Jain  1992 

Comprehensive  Reference 

Trescott  and  Larson  (1977),  Ames 
(1977),  Mitchell  and  Griffiths 
(1980),  Lapidus  and  Pinder  (1982) 

Computer  Program 

Trescott  et  al.  (1976),  Konikow  and 
Bredehoft  (1978),  Healy  (1990) 

Winter  (1983,  1986,  1989),  Sacks  (1992)  and  Krabbenhoft  and  Anderson  (1990) 
conducted  numerical  investigations  on  the  interactions  between  groundwaters  and  lakes. 
Both  Winter  and  Krabbenhoft  and  Anderson  used  fixed  finite-difference  meshes.  For  lakes, 
a  finite-difference  mesh  may  be  suitable;  however,  because  the  bed  of  most  wetlands  exhibits 
significant  slope,  a  poorly  designed  fixed  finite-difference  mesh  may  not  adequately 
characterize  water  fluxes  across  a  bed  situated  between  an  overlying  body  of  water  and  an 
underlying  aquifer.  The  issue  of  properly  modeling  fluxes  accross  a  sloping  bed  may  not  be 
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as  important  as  the  fact  that  accurate  simulations  of  wetland-flow  systems  require 
measurements  of  hydrologic  properties  at  finer  scales  and  frequencies  than  are  commonly 
made  (Hunt  1996,  Hunt  and  others  1996  1997). 

Assuming  the  hydrological  model  is  correct,  the  accuracy  of  any  computer  simulated 
hydrologic  system,  such  as  a  wetland,  depends  on  the  quality  of  the  data  used  to  calibrate 
model.  Through  a  well  instrumented  study-site  in  the  Kankakee  River  Valley  in  northwest 
Indiana,  the  USGS  and  the  USEPA  conducted  detailed  measurements  of  many  various 
variables  required  to  characterize  a  wetland  hydrologic  system.  McDonald  and  Harbaugh 
(1998)  developed  a  groundwater  model  using  the  available  data  and  conducted  the 
simulations  of  observed  hydrologic  process  using  the  computer  code  commonly  known  as 
MODFLOW.  The  focus  of  the  research  effort  was  to  study  the  effects  of  various  wetland 
restoration  activities  on  local  groundwater  levels.  Computer  simulations  were  used  to  predict 
the  effects  of  wetland  restoration  on  the  existing  hydrology  of  the  wetlands  and  adjacent 
farmlands  and  to  suggest  residence  times  and  path  ways  for  surface  water  and  groundwaters. 

The  finite-element  method  was  first  used  to  solve  groundwater  flow  and  solute 
transport  problems  in  the  early  1970s.  The  method  has  several  advantages:  (1)  irregular  or 
curved  aquifer  boundaries,  anisotropic  and  heterogeneous  aquifer  properties,  and  sloping  soil 
and  rock  layers  can  be  incorporated  easily  into  the  numerical  model;  and  (2)  the  accuracy  of 
solution  to  groundwater  flow  and  solute  transport  problems  is  quite  good.  The  main 
disadvantages  of  the  finite-element  method  include:  (1)  for  simple  problems,  the  finite 
element  method  requires  a  greater  amount  of  mathematical  and  computer  programming 
sophistication  than  does  the  finite-difference  method;  and  (2)  there  are  fewer  well- 
documented  computer  programs. 
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The  finite-element  method  provides  approximations  of  much  higher  order  than  does 
finite-difference  methodology.  The  whole  modeling  domain  can  be  subdivided  into  irregular 
elements  as  dictated  by  the  physical  geometry  of  the  problem.  The  size  of  each  element  can 
vary.  Small  elements  can  be  used  in  areas  where  rapid  change  occurs  with  model  parameters 
and  predicted  state  variables,  whereas  the  large  elements  can  be  used  where  such  variations 
are  less  severe.  Inhomogeneities  and  anisotropy  are  easily  accommodated  with  the  finite- 
element  method.  The  advantage  of  irregular  elements  and  higher-order  approximations  is 
certainly  important  when  the  "space"  domain  is  complex.  However,  this  advantage  is  less 
obvious  in  problems  involving  time  (evolution)  (Gupta  1975).  Selected  references  of  early 
finite-element  applications  to  model  groundwater  flow  are  listed  in  Table  3-2. 

3.3  Three-Dimensional  Deformable  Finite-Element 
Saturated  Groundwater  Model 

To  simulate  three-dimensional  groundwater  flow  in  phreatic  aquifer  systems  and  the 

exchange  of  water  between  a  wetland  and  the  contiguous  aquifer,  a  numerical  model  is 

needed  that  is  capable  of  modeling  the  free  groundwater  surface.  This  requires  first  that  the 

model  assume  a  tentative  location  of  the  free  surface  and  then  solve  for  the  distribution  of 

the  dependent  variable,  iJj,  over  a  fixed  domain.  All  this  requires  that  the  flux  boundary 

condition  across  the  surface  be  satisfied  and  that  the  pressure  be  atmospheric.  To  model  the 

exchange  of  water  between  a  wetland  and  an  aquifer,  it  is  necessary  to  recognize  that 

significant  vertical  and  lateral  movements  of  the  wetland  boundary  typically  occur  with 

changes  in  wetland  storage  (see  Figure  3-1).  This  transient  boundary  suggests  that  a 

deformable  finite-element  mesh  may  be  needed  in  the  model  to  characterize  adequately  the 

moving  groundwater/surface-water  interface. 
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Table  3-2.  Selected  references  of  early  finite-element  applications  to  model  groundwater. 


Topic 

Reference 

Early  Developments  in  Petroleum 
Reservoir  Modeling 

Price  etal.  (1968) 

Saturated  Groundwater  Flow 

Zienkiewicz  et  al.  (1966),  Javandel 
and  Witherspoon  (1968), 
Zienkiewicz  and  Parekh  (1970), 
Pinder  and  Frind  (1972). 

Unsaturated  Groundwater  Flow 

Neuman  (1973),  Gureghian  et  al. 
(1979),  Pickens  and  Gillham  (1980) 
Yeh  and  Ward  (1981),Huyakorn 
(1984),  Yeh  (1992). 

Solute  Transport 

Price  et  al.  (1968),  Guymon  et  al. 
(1970),  Neuman  (1973),  Van 
Genuchten  et  al.  (1977),  Kirkner  et 
al.  (1984). 

Application  to  Field  Problems 

Pinder  (1973),  Gupta  and  Tinji 
(1976),  Senger  and  Fogg  (1987). 

Comprehensive  References 

Ziekiewicz  (1971),  Pinder  and  Gray 
(1977),  Lapidus  and  Pinder  (1982), 
Huyakorn  and  Pinder  (1983). 

Computer  Programs 

Neuman  and  Witherspoon  (1970), 
Reeves  and  Duguid  (1975),  Segol  et 
al.  (1975),  Pickens  et  al.  (1979) 
Huyakorn  (1984),  Yeh  (1992). 

Numerical  groundwater  modeling  that  is  focused  on  simulating  the  interaction 
between  surface-water  and  groundwater  systems,  where  the  groundwater/surface-water 
interface  moves,  must  use  well-defined  linking  conditions  to  describe  the  exchange  of  water 
across  the  interface.  These  linking  conditions  are  unknown  in  the  beginning  of  the  modeling 
time  step  and  found  through  a  flux  type  interaction  that  then  forecasts  the  linking  boundary 
conditions  at  the  end  of  each  time  step.  Thus,  the  model  has  transient  boundary  conditions. 
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Wetland 


„  Water  Surface  at  Time  t 


Phreatic  Surface  at  Time  t,  Phreatic  Surface  at  Time  t. 

Groundwater  Aquifer 
Impermeable  Aquacludc 


Figure  3-1.  Conceptual  profile  model  of  a  wetland/groundwater  system 
Finally,  typical  finite-difference  modeling  is  not  able  to  trace  the  curvature  of  the  bottom  of 
a  wetland  (see  Figure  3-2).  Finite  elements,  however,  are  easily  configured  to  fit  non-linear 
spatial  geometries.  In  conclusion,  to  cope  with  the  above  described  three  technical  problems 
of  modeling,  i.e.,  the  free  surface  location,  incorporating  dynamic  changes  in  the  vertical  and 
horizontal  extent  of  the  ground/surface-water  interface,  and  tracing  the  curvature  of  a 
wetland  bottom,  a  special  three-dimensional  finite-element  model  is  needed. 


Finite-Element  Method 


Finite-Difference  Method 


Figure  3-2.  Comparison  of  finite-element  mesh  and  finite-difference  mesh 
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3.3.1  Governing  Equation  and  Boundary  Conditions 

Eq.  3-1  below  can  be  used  to  describe  incompressible,  saturated  groundwater  flow 
beneath  a  wetland  and  in  a  phreatic  aquifer: 


V{KVh)±Q  =  StQ    i,j  =  \,2,3 
dt 


(3-1) 


where  K  is  the  saturated  hydraulic  conductivity  tensor  (  LT  "');  h  is  equal  to  z+i|j  (L);  ijj  is 
pressure  head  (L);  z  is  elevation  head  ( L);  and  Ss  is  the  specific  storativity  (L*1).  From  Figure 
3-3,  the  boundary  of  volume  Q  is  characterized  through  a  combination  of  specified  flux 
boundaries  Tl  and  specified  head  boundaries  T2.  Pertinent  specific  no-flux  boundaries  include 
the  groundwater  phreatic  surface  (EC  and  DF)  and  the  bottom  impermeable  boundary  (GH). 
Relevant  specified  head  boundaries  T2  include  the  wetland  head  boundary  (AB  and  CD)  and 
lateral  boundaries  (EG  and  FH). 
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EG,  FH,  ACDB,  and  CD  are  specified  head  boundaries  T2 


Figure  3-3.  Setup  of  the  saturated  groundwater  problem  around  a  wetland 
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As  to  the  boundary  condition  controlled  by  the  wetland  water  table,  the  initial  water 
table  of  the  wetland  and  the  pumping  discharge  were  taken  as  the  only  inputs  of  the  coupling 
model.  During  pumping,  the  pumping  discharge  lowered  the  wetland  water  table,  and  it  also 
triggered  an  inflow  from  the  aquifer  to  the  wetland  through  the  peat  layer.  In  each  simulation 
time  step,  the  pumping  discharge,  estimated  inflow  from  the  aquifer  to  the  wetland,  initial 
water  table  of  the  wetland,  and  relation  between  the  volume  of  the  wetland  and  wetland 
water  table  were  used  as  inputs  in  a  mass  balance  equation  that  was  solved  to  determine  the 
water  table  at  the  end  of  the  time  step.  The  estimated  water  table  at  the  end  of  the  time  step 
was  put  back  into  the  saturated  groundwater  model  to  adjust  the  inflow  from  the  aquifer  to 
the  wetland  estimated  initially.  Therefore,  several  iterations  in  each  time  step  were  necessary 
to  get  an  accurate  estimated  inflow  from  the  aquifer  to  the  wetland  and  a  convergent  solution 
of  the  water  table  in  the  wetland  at  the  end  of  the  time  step.  When  pumping  was  stopped,  the 
pumping  discharge  was  taken  as  zero.  Figure  3-4  shows  the  procedure  of  finding  the  inflow 
from  the  aquifer  to  the  wetland  and  the  coupling  boundary. 
3.3.2  Finite  Element  Model  of  Saturated  Groundwater  Flow 

To  obtain  a  finite  element  approximation  of  Eq.  3- 1 ,  two  spaces,  5  and  V,  are  defined 
(Hughes  1987): 

S=  {u\ueH\u(xyj)\T=qr^  (32) 

V=  {u\ueHl,u(x,y,z)\r^o}  (3.3) 

where  H1  is  Sobolev  space  comprised  of  functions  which  are  square-integrable  generalized 
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Wetland 
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Simulate  groundwater  flow  based 
on  the  predicted  wetland  water  table 
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Find  the  groundwater  inflow  to  the  wetland 
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Figure  3-4.  Flowchart  illustrating  the  iteration  procedure  of  finding 
the  coupling  boundary  condition 

derivatives  through  order  1;  u  is  a  space  function  that  is  member  of  Sobolev  space;  qn  is  the 

vertical  flux  along  r2.  and  x,  y,  z  are  three  spatial  coordinates.  Any  function  from  S  space  is 

equal  to  the  flux  boundary  value  along  specified  head  boundary  (r2).  Also,  any  function  u 

from  V  space  is  equal  to  zero  along  head  boundary  (r2).  Assuming  a  basis  function  w  is 

chosen  from  V  space,  then  the  following  integration  of  Eq.  3-1  over  space  domain  Q  is 

always  true: 
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w[V-{KVh)]dQ±\    wQdQ=l  wSs—dQ,    i,j=l,2,3  (3.4) 


where  Eq.  3-4  expresses  the  governing  groundwater  flow  equation  integrated  over  the 
modeling  domain. 

Manipulation  of  the  first  term  on  the  left  hand  side  of  Eq.  3-4  begins  with  the 
application  of  Green's  theorem,  which  is  an  expanded  form  of  the  divergence  theorem.  This 
converts  the  integral  into  two  terms,  one  of  which  is  evaluated  only  on  the  surface  of  the 
region  to  be  simulated.  Green's  theorem  is: 


(V-P)AdQ  =  l  (Pn)AdT-  I  (P-VA)dQ 

J  J  0-5) 

a  r  a 


where  A  is  a  scalar,  and  P  is  a  vector  quantity.  The  boundary  of  volume  Q  is  defined  by  T, 
and  n  is  a  unit  outward  normal  vector  to  the  boundary.  Application  of  the  Galerkin  criterion 
(Hughes  1986)  and  Green's  theorem  to  Eq.  3-4  leads  to: 


w[(KVh)-n]dT+\    w[{KV h)-n]dT- \  (Vw)-{KVh)dQ 

Jr2  «/Q 


±1  w-QdQ=l   wS  —  dCl 

\a      Sdt  (3-6) 


s 


31 

Because  w  is  a  member  of  V  space,  it  is  equal  to  zero  at  all  specified  head  boundaries  (r2); 
thus,  the  second  term  in  Eq.  3-6  equals  zero.  If  it  is  further  assumed  that  Q  equals  zero,  then 
Eq.  3-6  reduces  to: 


(    w[{KVh)-n]dT-\  {Vw)-(KVh)dQ=  I   wSg — -dQ 
Jr,  t/Q  Jq 


(3-7) 


For  an  approximate  solution  of  h(x,  y,  z,  t),  it  is  assumed  that: 


^W(x^,z,0^^(eW,2)^(0  ,  NteV  ,    heS  (3-8) 


/  =  i 


■t^e\x,y,: 


ww=Z,N?Xx,y,z)ct    ,  NteV  (3-9) 

;=1 


where  h(e)  is  the  approximate  solution  for  hydraulic  head  within  element  e,  Nj(e)  are 
interpolation  functions  for  each  node  within  element  e,  p  is  the  number  of  nodes  within 
element  e,  and  d;(t)  are  the  unknown  values  of  hydraulic  head  for  each  node  within  element 
e  and  at  time  t.  At  the  specified  head  boundaries,  T2,  dj(t)  =  h(x,y,z)|p2  Wj<e)  is  the  weighting 
function  for  node  i,  and  the  limits  of  the  integration  are  chosen  to  represent  the  volume  of  the 
elements.  In  Galerkin's  method,  the  weighting  function  for  each  node  in  the  element  is  set 
equal  to  the  interpolation  function  for  the  node,  Wj(e)  =  Nj(e). 

When  the  approximate  solutions  (Eq.  3-8  and  Eq.  3-9)  are  substituted  into  Eq.  3-7, 
the  resulting  equation  is  not  satisfied  exactly;  consequently,  an  error  or  residual  occurs  at 
every  point  in  the  problem  domain.  However,  if  it  is  also  assumed  that  the  aquifer  is 
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isotropic,  that  is  to  say,  values  of  saturated  hydraulic  conductivity  in  the  three  coordinate 
directions  are  constant  within  an  element  (but  can  vary  from  one  element  to  the  next),  the 
contribution  of  any  element  e  to  the  residual  at  a  node  i  associated  with  the  element  can  be 
described  through  the  following  system  of  equations: 


R}e\t) 


R':\t) 
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rW 


dy® 
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(3-10) 


where  p  =  the  number  of  nodes  in  each  element  (i.e.,  triangle  element:  p  =  3),  and  Rj(e)  =  the 
residual  error  at  node  i  in  element  e: 
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where  V(e)  is  the  volume  of  element  e.  Eq.  3-10  represents  a  system  of  equations  pertinent 
to  a  single  element,  and  there  is  an  equation  for  every  node  in  the  element.  When  equations 
for  all  elements  in  the  mesh  are  combined,  the  global  equations  are  obtained  for  the  system: 


*i(0 

\v\ 

rf,(0 

*»(') 

-  [K\ 

4,(0 

(3-14) 


By  setting  the  residuals  equal  to  zero,  we  have: 
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*\C 


ddAt) 


dt 


dd„(t) 


dt 


[F] 


(3-15) 


Eq.  3-15  is  a  system  of  ordinary  differential  equations,  whose  solution  provides  values  of 
head  d  and  the  derivative  with  respect  to  time,  dd/dt,  at  each  node  in  the  finite-element  mesh. 
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Although  several  methods  are  available  for  solving  this  system  of  equations,  it  has  become 
standard  practice  to  use  the  finite-difference  method: 


dh,  ,      d(t+At)-d(t) 
-(e)  =  -^ — /     w  (3-16) 


dt  At 


If  a  variable,  oo  is  defined  ,  such  that: 


t-t 

G>  = 

At 


Then: 


(3-17) 


d(e)=  (l-u)J(0+o)^+AO  (3-18) 

which  can  be  extended  to  the  vector  of  [d]: 

[d(e)]  =  (1  -Q^Oj+w  [d(t+At)}  (3-19) 

If  Eq.  3-19  and  Eq.  3-18  are  substituted  into  Eq.  3-15,  the  finite-difference  formulation  for 
the  transient  saturated  flow  equation  is  as  follows: 

;  [C]  +  co  At[K]  )[d(t+At)}  =(  [C]  -  ( 1  -w)  At[K\) [d(t)]  +A/((1  -a>)[F(t)]+vJ[F(t+At)] 


(3-20) 
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Eq.  3-20  represents  a  system  of  linear  equations  with  unknown  variables  of  d;(t),  and 
where  Nj(x,y,z)  are  the  basis  functions  in  a  tri-linear  hexahedral  element  domain.  Nj(x,y,z) 
does  not  change  with  time.  The  most  common  basis  functions  used  are  hat  functions,  which 
are  linear  functions  of  space  variables. 

The  solution  of  Eq.  3-20  begins  by  specifying  the  initial  values  of  [d]  (i.e.,  the  values 
of  head  at  time  t0=  0).  Next,  the  system  Eq.  3-20  is  solved  to  obtain  values  of  {d(t+At)}  at 
the  end  of  the  first  time  step  1.  This  solution  process  is  repeated  for  the  next  time  step,  and 
so  on  until  the  specified  time  duration  of  simulation  is  satisfied.  Depending  on  the  choice  of 
co,  several  different  subsets  of  the  finite-difference  formulations  are  defined. 
When  co=0,  it  is  the  forward-difference  method  is  obtained: 

(  [C]  )[d(t+At)}  =(  [C]  -  At[K\)[d(t)}  +Ar[F(0]  (3-21) 

When  co=l/2,  the  central-difference,  or  Crank-Nicholson,  method  is  obtained: 

[C]  +  1a*M  WAf)]=(  [C]-(  l-±)to[K\)[d(t)]  +  ¥{F(t)]+[F(t)]) 


2 


2) 


2 


(3-22) 


When  co=l,  it  is  the  backward-difference  method  is  obatined: 


(  [C]  +  At[K]  ){d(t+At)}  =  [C]  [d(t)]  +A/[F(0] 


(3-23) 
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3.3.3  Adaptation  of  the  Saturated  Finite-Element  Model  to  Describe  the  Free  Surface 

Several  papers  have  been  written  on  the  groundwater  free-surface  problem 
(e.g.,Taylor  1967,  Neuman  1970,  Finn  1976,  France  1976,  and  Huyakorn  1986).  Most  of 
these  papers  discussed  the  steady-state  free-surface  problem  except  for  France  (1976).  This 
thesis  examines  the  transient  problem  starting  from  the  basic  ideas  that  Neuman  (1970) 
presented  for  steady  flow.  Basically ,  both  France  (1976)  and  Neuman  (1970)  tried  to  satisfy 
the  two  boundary  conditions  and  are  consistent.  This  chapter  is  also  based  on  the  two 
conditions  and  adapts  the  method  they  used  to  locate  the  free  surface. 

The  primary  dependent  variable  to  be  used  in  modeling  the  free  surface  presented 
here  is  the  hydraulic  head  h  defined  as: 


h  =  z  +  ^  (3-24) 

Y 


where  z  is  the  vertical  coordinate;  y  is  specific  weight  of  a  unit  volume  of  water  pg  and  v|/ 
is  the  pressure.  Along  the  free  surface,  there  are  two  conditions  which  must  be  satisfied. 
First,  i|/=0  ( pressure  atmospheric),  and  hence  for  a  given  position,  h  is  prescribed  on  it.  That 
is: 

k=Z  (3-25) 

For  the  steady-state  scenario,  the  second  condition  is  simply  that  there  is  a  specified 
flux  across  the  free  surface  boundary  (i.e.,  a  zero  flux  or  a  flow  imposed  equal  to  net 
recharge).  Thus,  in  addition  to  Eq.  3-24  at  the  free  surface,  there  is: 
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K®L-i+r§hj+K§!Lk  =  q  (3-26) 

xdx         ydyJ        zdz        Ho 


If  steady  state  has  not  been  reached  and  the  boundary  is  moving  with  a  velocity 
normal  to  its  instantaneous  configuration  of  Vn,  then  the  quantity  of  liquid  entering  its  unit 
area  is  modified,  and: 


dh  -.     „  dh  -.     „  dh  - 
—  i  +  K  —  j  +  K  — 

dx         ydyJ       zdz 


K*—.  l  +  Ky^7.J  +K^k  =  ^o  +  VnS  -  %  (3-27) 


where  i,  j,  k  are  unit  orthogonal  vectors;  S  is  the  specific  yield  coefficient  (or  the  effective 
porosity)  relating  the  total  volume  of  material  to  the  quantity  of  liquid  which  can  be  drained 
from  it,  and  q^  equals  the  net  flux  across  the  free  surface. 

The  transient  problem  is  solved  by  considering  the  solution  to  be  a  number  of  steady- 
state  problems  at  small  intervals  of  time  At.  For  a  given  time  level,  the  governing  equation 
is  solved  to  obtain  values  of  nodal  heads  and  velocities  of  free  surface.  The  components  of 
actual  velocities  are  determined  by  dividing  by  the  specific  yield  coefficient.  By  repeating 
this  process  for  all  modal  points  on  the  free  surface,  a  new  configuration  and  a  new  set  of 
boundary  conditions  are  defined.  The  actual  nodal  coordinates  of  the  free  surface  are  found 
using  the  Newton-Raphson  iterative  method.  Thus,  a  new  finite-element  mesh  is  generated 
by  modifying  the  coordinates  of  the  nodes  in  the  new  flow  domain. 
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3.4  Validating  the  Code  of  Three-Dimensional  Deformable 
Finite-Element  Groundwater  Flow  Model 

3.4.1  Theis  Solution 

To  validate  the  code  of  the  three-dimensional  deformable  finite-element  groundwater 
flow  model,  a  comparison  was  made  of  results  generated  by  the  numerical  model  and  those 
created  with  an  analytical  model  describing  unsteady  flow  to  a  well  in  a  phreatic  aquifer. 
Figure  3-5  shows  a  schematic  of  unsteady  flow  to  a  well  in  a  phreatic  aquifer. 
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Figure  3-5.  Unsteady  flow  to  a  well  in  a  phreatic  aquifer  at  time  t. 
where  s  =  drawdown(L),  H0=  initial  saturated  aquifer  thickness  (L),  Q  =  pumping  discharge 
(L3/T),  and  r  =  radial  distance  from  the  well  (L). 

As  shown  in  Figure  3-2,  the  boundary  condition  along  the  phreatic  surface  is 
nonlinear  and  the  shape  and  position  of  the  phreatic  surface  are  unknown.  As  a  result,  an 
exact  analytical  solution  that  describes  groundwater  flow  is  not  known.  However,  by  the 
employing  the  Dupuit  assumptions  and  assuming  recharge  is  zero,  the  unsteady  flow  to  a 
fully  penetrating  well  in  a  phreatic  aquifer  can  be  described  in  radial  coordinates  by: 


d2h2  +\M2  _  2S  dh 
dr2     r    dr        K   dt 


(3-28) 
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Initial  condition:  h(t0)  =  H0 

Boundary  condition:  h(t)  |  {I=m)  =  H0 

Here,  S  is  the  phreatic  aquifer  specific  yield  (or  effective  porosity,  ne).  If  drawdowns  are 
sufficiently  small  compared  to  the  average  depth  of  flow,  (say  0.02Ho),  then  the  approximate 
solution  to  this  problem  is  the  Theis  solution  (Bear  1976): 

4%T  (3-29) 

where  u  =  SS/(4T  t);  W(u)  is  the  well  function,  and  T=KH0=  initial  transmissivity. 

The  conditions  that  define  the  test  problem  used  to  validate  the  numerical 
groundwater  flow  model  include:  a  pumping  discharge,  Q,  of  1.5  m3/hr;  an  initial  saturated 
phreatic  aquifer  thickness,  H0,  of  20  m;  a  specific  yield  S  (or  ne),  of  0.1;  and  a  hydraulic 
conductivity,  K,  of  15  m/day.  For  the  three-dimensional  finite -element  model,  the  specific 
storativity,  Ss,  is  0.0001  m"1;  where  Ss=pg(a+nP);  n  equals  porosity;  a  equals  solid  matrix 
compressibility  [LT"2M-1];  g  is  acceleration  of  gravity  [LT"2];  and  P  is  fluid  compressibility 
[LT2M"'].  Because  of  problem  symetry,  modeling  was  limited  to  simulating  hydraulic  head 
variations  over  a  wedge  section  of  the  aquifer  using  a  21x1 1x3  mesh  (Figure  3-6). 

The  boundary  condition  controlled  by  well  water  table,  the  initial  water  table  of  well 
and  the  pumping  discharge  were  taken  as  the  only  input.  The  pumping  pulled  water  out  of 
the  well,  on  the  other  hand,  the  groundwater  supplied  well.  In  each  simulation  time  step, 
based  on  the  known  water  table  in  the  beginning  of  time  step,  and  an  estimated  water  table 
of  well,  a  new  discharge  supplied  from  groundwater  to  well  was  calculated  from  the 
groundwater  modeling.  Based  this  new  calculated  supplying  discharge  to  the  well,  the 
pumping  discharge,  and  the  circular  area  of  the  well,  a  more  accurate  well  water  table  could 
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be  estimated  again.  The  groundwater  flow  model  was  solved  again.  After  many  iterations, 
a  convergent  and  accurate  well  water  table  could  be  obtained  before  the  modeling  of  next 
time  step. 
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Figure  3-6.  Finite-element  mesh  of  the  wedge  used  to  simulate  of 
a  well  pumping  in  phreatic  aquifer 

Because  the  mesh  is  deformable,  the  known  pressure  at  the  free  surface  and  flux  were 
used  to  locate  the  free  surface.  During  calculation,  a  slice  successive  over-relaxation  (SSOR) 
method  was  used  as  an  iterative  approach  to  locate  the  free  surface.  The  results  illustrated 
in  Figure  3-7  validate  the  coding  of  the  numerical  model. 
3.4.2  Groundwater  Flow  around  a  Circular  Wetland 

To  validate  the  three-dimensional  finite-element  saturated  groundwater  flow  model 
again,  the  numerical  model  was  used  to  simulate  the  groundwater  flow  around  a  circular 
wetland  and  induced  by  a  flood  wave  in  the  wetland.  It  is  assumed  that  the  wetland  stage 
fluctuations  can  be  described  by  a  nonperiodic  Gaussian  fucntion  of  time  as  Eq.  3-30: 
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Figure  3-7.  Comparison  of  the  numerical  and  analytically  simulated  head 
for  the  example  problem  of  a  transient  pump  test  in  a  phreatic  aquifer 


Unlike  the  coupling  boundary  condition  that  was  used  in  the  last  case,  a  flood 

wave  boundary  condition  used  in  this  numerical  model  here.  The  solid  line  and  the  dash 

line  in  Figure  3-8  are  the  three-dimensional  finite-element  model  simulated  results.  The 

triangle  points  and  the  diamond  points  are  the  groundwater  analytical  solutions  from  Eq. 

2-24  obtained  in  Chapter  2.  The  circular  wetland  which  has  a  radius  of  50  meter  was 

assumed.  All  other  system  parameters  are  given  in  the  Table  3-3.  The  numerical  model 

and  analytical  solution  matched  very  well. 

3.5  A  Numerical  Model  Application  to  Simulate  the  Interaction 
Between  a  Wetland  and  an  Aquifer 

The  finite-element  model  was  applied  to  simulate  the  interaction  between  a  wetland 

and  a  contiguous  phreatic  aquifer  stressed  by  an  active  productive  well.  A  horizontal, 

irregular  quadrilateral  mesh  was  used  to  define  the  irregular  boundaries  of  the  wetland  and 

the  area  beyond  the  wetland.  Figure  3-9  illustrates  the  numerical  mesh  corresponding  to  the 
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Figure  3-8.  Comparison  of  the  numerically  and  analytically  simulated  groundwater  table. 

Table  3-3.  Parameters  and  their  values 


Parameter 

Value 

Parameter 

Value 

H00  (m) 

100.0 

T    (m2/h) 

0.635 

Ar   (m) 

0.25 

S 

0.01 

At    (hr) 

0.02 

t0     (h) 

5.0 

k    (m/hr) 

0.00635 

o 

0.4 

horizontal  and  vertical  extent  of  the  phreatic  aquifer.  Beyond  the  wetland,  the  vertical  extent 
of  the  mesh  defines  the  thickness  of  the  phreatic  aquifer.  Over  the  inundated  portions  of 
wetland,  the  thickness  of  the  mesh  defines  the  vertical  distance  between  the  wetland  bottom 
and  the  underlying  aquiclude.  The  wetland  is  shown  in  the  middle  of  this  area.  The  model 
domain  is  104x104x10  meters.  The  four  lateral  boundaries  and  bottom  boundary  are 
assumed  to  be  impermeable.  An  active  production  well  is  located  within  the  model  domain. 
It  was  assumed  that  the  initial  groundwater  table  and  the  water  stage  of  wetland  were  the 
same.  The  deepest  water  depth  in  the  wetland  is  1.0  meter.  Model  runs  were  conducted  to 
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simulate  the  exchange  of  water  between  wetland  and  aquifer.  Such  interactions  between  the 
wetland  and  the  aquifer  were  examined  under  transient  pumping  conditions  and  precipitation. 
Since  the  real-time  stage  of  the  wetland  is  used  as  a  boundary  condition  over  the  wetted  part 
of  wetland,  the  horizontal  and  vertical  extent  of  the  mesh  must  be  adjusted  to  trace  the 
vertical  and  horizontal  movement  of  the  transient  wetland  water  surface. 
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Upper-extent  of  the  phreatic  aquifer  defined 
the  bottom  elevation  inundated  wetland 


Pumping  well 


aquiclude 


5-2005E-04 
1  9999E+01 
3.9999E+01 
59999E+OI 
7.9999E+01 
9  9999E+0I 


Figure  3-9.  Three-dimensional  finite-element  mesh  for  the  simulation  of  wetland  and  aquifer 
The  modeled  area  was  divided  into  a  21x21x10  hexahedron  mesh.  In  the  middle,  a 
9x9  quadrilateral  mesh  was  used  to  represent  the  wetland  area  (Figure  3-9).  Hydrologic 
conditions  that  define  the  test  problem  include:  a  net  precipitation  rate  of  0.02  m/hour;  a 
specific  yield,  S,  of  0.075.  a  hydraulic  conductivity,  K,  of  10  m/day;  and  a  pumping 
discharge,  Q,  of  1.554  nvVhr.  Finally,  data  from  Figures  3-10  and  Figure  3-11  were  used  to 
define  the  relationship  between  the  wetland  water  stage,  the  surface  water  area  of  wetland, 
and  the  corresponding  wetland  volume. 

Figure  3-12  illustrates  the  initial  condition  where  it  was  assumed  that  the  wetland 
surface-water  stage  and  phreatic  surface  elevation  were  the  same.  Figure  3-13  illustrates  a 
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plan  view  of  the  numerical  mesh  at  time  zero.  Figure  3-14  through  Figure  3-21  represent 
various  depictions  of  head  contours  and  iso-surface  plots  of  head  predicted  at  hour  5  of  the 
simulation.  These  figures  show  that  water  moves  laterally  into  the  wetland  and  pumping  well 
because  the  same  net  rainfall  forces  the  water  table  to  increase  faster  than  in  the  water  stage 
in  the  wetland.  The  directions  of  flow  also  depend  on  the  active  production  well.  If 
production  is  large,  more  water  flows  to  the  well  instead  of  flowing  to  the  wetland. 
Otherwise,  some  water  flows  towards  and  into  the  wetland  as  groundwater  inflow.  Beneath 
the  wetland,  the  hydraulic  head  increases  with  time  and  with  increasing  depth  (Figure  3-17 
and  Figure  3-18),  which  shows  that  water  flows  from  the  aquifer  into  the  wetland.  As  to  the 
pumping  well,  the  contours  are  almost  vertical  and  parallel  to  each  other  near  the  well 
(Figure  3-20  and  Figure  3-21).  Groundwater  moves  towards  the  well  in  a  horizontal  direction 
because  the  well  fully  penetrates  the  aquifer. 
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Figure  3-10.  Relationship  between  the  area  of  wetland  water  surface. 
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Figure  3-11.  Relationship  between  the  deepest  water  depth  and  the  volume  of  wetland. 


Figure  3-12.  Dlustration  of  the  initial  condition  where  it  was  assumed  the  wetland  surface 
water  stage  and  the  phreatic  surface  elevation  were  the  same 
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Figure  3-13.  Plan  view  of  the  numerical  mesh  at  time  zero. 
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Figure  3-14.  Three-dimensional  finite-element  mesh 
and  simulated  head  contours  at  hour  5. 
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Figure  3-15.  simulated  head  contour  at  hour  5. 


Figure  3-16.  Simulated  iso-surface  plot  of  10.6  meter  head  at  hour  5. 
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Figure  3-17.  Three-dimensional  finite-element  mesh  and  head  contour 
on  a  cross  section  through  wetland 
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Figure  3-18.  Simulated  head  contour  on  the  vertical  cross  section  through  the  wetland 
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Figure  3-19.  Modeled  iso-surface  plot  of  head  equal  to  1 1.1  m,  1 1  m  and  10.5  m 
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Figure  3-20.  Mesh  and  modeled  head  contour  with  vertical  cross  section 
through  the  pumping  well 
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Figure  3-21.  Head  contour  (m)  on  the  vertical  cross  section  through  the  pumping  well 

A  computer  with  a  Pentium  1.4GHz  processor  took  4  hours  to  complete  the  5-hour 

simulation  using  a  time  step  of  0.1  hours  and  3  iterations  per  time  step.  To  use  this  model 

over  a  real  system,  a  finer  mesh  would  be  desirable.  Unlike  most  rectangular  finite-difference 

grids,  the  finite-element  grid  used  here  precludes  the  use  of  more  efficient  solution 

algorithms  that  divide  simulation  domain  into  three  individual  directions.  In  order  to  use  this 

model  to  investigate  wetland  and  aquifer  interaction  on  a  large  scale,  it  would  be  desirable 

to  introduce  a  simplified  mesh  that  would  in  turn  be  more  amenable  to  efficient  solvers.  With 

computer  hardware  developing,  the  simulation  speed  can  be  accelerated  consequently.  To 

improve  the  simulation  speed  of  the  model  from  fundamental,  a  kind  of  hybrid  method  which 

combines  finite -element  method  and  finite-difference  method  could  be  an  effective  choice 

for  a  large  scale  practical  problem. 

3.6  Application  of  Three-Dimensional  Deformable  Saturated 
Groundwater  Flow  Model  to  SV5  Wetland 

3.6.1  Site  Description-SV5 

SV5  is  an  isolated  wetland  found  in  pine  flatwoods  located  in  southeast  Florida. 
It  is  part  of  the  Savannas  wetland  system  located  south  of  the  city  of  Jensen  Beach  in  Martin 
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County.  SV5  is  located  approximately  two  km  west  of  the  Atlantic  Ocean  and  is  adjacent  to 
a  large  municipal  well  field.  The  South  Florida  Water  Management  District  (SFWMD) 
monitors  this  wetland  system.  There  are  13  surficial  aquifer  municipal  water-supply  wells 
located  from  a  few  hundred  meters  to  four  kilometers  from  SV5  (Figure  3-22).  The  wetland 
is  roughly  circular  in  shape,  with  a  diameter  of  approximately  60  meters  (Wise  et  al.  2000, 
Walser  1998). 
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Figure  3-22.  Well  Locations  (Switt  et  al.  1998) 
The  surficial  aquifer  that  underlies  the  wetland  is  composed  of  three  layers.  The  first 
layer,  which  is  located  from  land  surface  to  12- 18m  below  the  surface,  consists  of  white,  gray 
and  brown,  largely  fine-  to  coarse-grained  quartz  sand  intermixed  with  shell  beds.  It  receives 
recharge  directly  from  precipitation.  The  second  layer,  which  is  3-6  m  thick,  is  comprised  of 
fine  to  very  fine  sand  with  small  amounts  of  shell  and  clay.  The  third  layer,  which  extends 
to  a  depth  of  between  37-43  m,  consists  of  limestone  intermixed  with  sand  and  shell.  The 
upper  confining  layer  lies  beneath  the  surficial  aquifer.  This  layer  consists  of  low- 
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permeability  rocks  that  are  primarily  clastic.  The  Floridian  aquifer  system  lies  under  the 
upper  confining  layer.  The  system  is  quite  thick  and  consists  of  limestone  of  variable 
permeability  (Fetter  1994).  The  climate  in  the  area  is  considered  humid  sub-tropical,  and  the 
area  receives  an  average  of  1.5  m  of  rainfall  per  year  (Walser  1998). 
3.6.2  Field  Methods  and  Wetland  Aquifer  Interaction  Test 

The  SV5  wetland  aquifer  interaction  test  was  conducted  (Wise  et  al.  2000).  Eighteen 
monitoring  wells  were  installed.  Wells  were  placed  along  transects  running  from  east  to  west 
and  north  to  south.  Six  wells  were  placed  outside  the  wetland  on  the  western  transect  and  two 
on  each  of  the  other  three  transects.  The  wells  were  identified  by  the  transect  designation 
followed  by  the  distance  location  from  the  center  of  the  wetland,  i.e.,  well  37  is  the  well 
located  37  m  from  the  center  of  the  wetland  along  the  west  transect.  SFWMD  installed  two 
wells  and  a  staff  gage  at  the  site.  The  wells  measure  groundwater.  The  groundwater  well  is 
screened  at  a  depth  of  six  meters  and  grouted  from  the  surface  to  a  depth  of  five  meters. 

The  water  and  peat  depths  in  the  wetland  were  measured  along  four  transects  through 
the  center  of  the  wetland  at  45°  angles  from  one  another.  Peat  and  water  depth  measurements 
were  taken  at  three-meter  increments.  The  peat  depth  was  determined  by  pushing  a  piece  of 
1  cm  rebar  into  the  peat  at  a  steady  rate  until  significant  resistance  was  felt.  This  method  is 
most  easily  applied  to  sites  such  as  SV5  that  are  free  of  significant  tree  roots.  The  depth  of 
the  water  table  was  determined  by  measuring  the  water  level  in  the  exterior  monitoring  wells. 
Elevation  data  for  the  sites  were  measured  utilizing  a  Topcon  Self-Leveling  Rotating  laser 
level.  The  South  Florida  Water  Management  District  staff  gage  on  the  site  was  used  as  a 
relative  reference  for  the  elevation  data,  and  all  elevation  data  were  expressed  in  terms  of 
NGVD  of  1929. 
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The  water  and  peat  depth  information  was  entered  into  Surfer,  a  three-dimensional 
graphical  interpolation  program,  to  develop  three-dimensional  contours  of  the  wetland.  A 
stage  -volume  curve  to  describe  the  wetland  was  developed  by  utilizing  Surfer  to  determine 
the  wetland  volume  at  varying  water  levels.  A  regression  curve  was  fitted  to  this  data  to 
determine  the  equation  that  describes  the  stage-volume  curve  (Switt  et  al.  1998,  Wise  et  al. 
2000).  This  curve  describes  the  top  of  the  peat  layer.  The  process  was  repeated  to  develop 
a  curve  that  describes  the  bottom  of  the  peat  layer. 

The  goal  of  the  test  was  to  stress  hydraulically  a  wetland  by  pumping  the  standing 
water  to  lower  the  water  level  a  predetermined  amount,  which  varies  from  wetland  to 
wetland,  and  then  monitoring  the  response  in  wells  and  subsequent  recovery  of  the  water  in 
the  wetland.  For  the  wetland  test  to  be  effective,  the  wetland  must  be  a  discharge  wetland  at 
the  end  of  the  pumping  phase.  The  wetland  test  was  conducted  utilizing  a  10.2  cm  centrifugal 
trash  pump  to  pump  the  water  from  the  wetland  until  the  water  level  in  the  wetland  dropped 
30.5  cm.  The  water  was  pumped  into  a  drainage  ditch  located  approximately  100  m  from  the 
wetland.  This  distance  was  chosen  to  ensure  that  the  water  pumped  from  the  wetland  did  not 
reenter  the  system.  The  discharge  rate  of  the  pump  was  approximately  70  m3/hr.  The 
discharge  rate  was  measured  by  determining  the  amount  of  time  it  took  to  fill  a  five-gallon 
plastic  bucket  at  the  end  of  the  discharge  hose.  It  took  nine  hours  to  lower  the  wetland  water 
level  30.5  cm. 

During  the  testing,  water  levels  in  all  of  wells  were  manually  measured  every  hour, 
with  the  exception  of  the  surface  water  well,  which  was  measured  every  30  minutes,  utilizing 
a  Slope  Indicator  water  tape.  Once  the  pump  was  turned  off,  water  levels  in  all  of  the  wells 
were  manually  measured  every  hour  for  22  hours.  Exterior  and  interior  well  water  levels 
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were  recorded  for  the  wells  in  the  wetland  area. 
3.6.3  Stage- Volume  Relationships 

Water  and  peat  depth  data  were  used  to  produce  a  three-dimensional  contour  plot  of 
SV5  (Figure  3-23  and  Figure  3-24).  The  regression  equation  that  describes  the  volume  of 
water  above  the  peat  layer  is 


V  =3454S2- 148.455 


(3-30) 


where  Vw  is  the  wetland  volume  (m3)  at  s,  the  stage  level  (m)  (Switt  et  al.  1998,  Wise  et  al. 
2000).  The  bottom  of  the  peat  layer  is  described  by  the  follow  regression  equation  (Figure 
3.25): 


V  =-116s4+607s3+938,s2  +  108s 
p 


(3,31) 


where  Vp  is  the  volume  above  the  bottom  of  the  peat  layer  (m3);  and  Vp-Vw  ~  volume  of 
saturated  peat.  The  coefficient  of  determination  equals  unity  for  both  of  these  regressions. 


Figure  3-23.  Top  surface  of  peat  layer  of  wetland 
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Figure  3-24.  Bottom  surface  of  peat  layer 
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Figure  3-25.  Relationship  between  wetland  volume  and  water  stage 


56 

Switt  et  al.  (1998)  used  a  similar  formula  (Eq.  3-32)  to  calculate  the  interaction 
between  a  wetland  and  the  underlying  groundwater: 


SP=K\  ^-^\SAW 

\     L    )      w  (3-32) 


where  SP=  exchanging  flux  (mVday),  Ks=  hydraulic  conductivity  (m/day),  G  =  groundwater 
level  (m),  W  =  wetland  water  level  (m),  L  =  average  peat  thickness  (m),  and  SA  w  =  wetland 
wetted  surface  area  (m2).  Over  the  de-watered  area  of  the  wetland,  seepage  decreases  from 
the  above  value  at  the  edge  of  the  water  to  zero  at  the  edge  of  the  wetland.  Switt  et  al.  (1998) 
assumed  a  standing  water  head  difference  equal  to  1/3  of  the  value  use  in  Equation  3-8.  Thus, 
seepage  for  the  de-watered  area  of  the  wetland  was 


SP=K 


(  G-W" 


SAm  (3-33) 


where  SANW=  wetland  dry  surface  area  (m2). 

The  hydraulic  conductivity  of  the  organic  soil  in  the  peat  layer  is  critical  in  wetland 
systems  because  it  controls  the  amount  of  groundwater  that  flows  from  the  aquifer  into  the 
wetland  or  from  the  wetland  into  the  aquifer  (Truesdell  1995).  Determining  a  value  of 
hydraulic  conductivity  for  the  soils  in  wetland  system  is  particularly  challenging. 
3.6.4  Creation  of  Finite-Element  Mesh 

The  model  domain  was  approximately  three  times  the  size  of  the  wetland.  Because 
the  wetland  is  approximately  radically  symmetrical  it  was  possible  to  develop  a  highly 
refined  4.0°  numerical  wedge  to   represent  the  wetland  aquifer  system.  The  aquifer  was 


57 
divided  into  10  layers  to  obtain  a  1 1x26x3  mesh  (Figure  3-26).  The  peat  layer  was  divided 
into  three  layers  producing  a  3x1 1x3  mesh.  The  mesh  was  divided  and  fitted  along  the  shape 
of  peat  layer. 


Wetland  Peat  Layer 


Figure  3-26.  Finite-element  mesh  of  wedge  of  wetland 
and  aquifer  for  the  saturated  flow  model 


3.6.5  Solution  Procedure 

To  simulate  the  wetland-aquifer  interaction  observed  at  SV5,  the  model  must 
correctly  simulate  the  movement  of  the  free  surface  as  changes  occur  with  the  wetland  water 
stage.  Two  boundary  conditions  must  be  satisfied  on  free  surfaces;  first,  there  can  be  the  flux 
across  the  surface  must  be  represented,  and,  second,  the  pressure  must  be  atmospheric.  The 
initial  fixed  domain  solution  usually  will  not  satisfy  both  boundary  conditions 
simultaneously.  Therefore,  the  trial  free  surface  is  adjusted  and  a  solution  obtained  for  the 
revised  fixed  domain.  This  process  is  repeated  until  the  required  adjustment  in  the  free 
surface  is  considered  negligible.  Thus,  the  solution  of  the  non-linear  free  surface  problem  is 
achieved  by  a  series  of  solutions  to  linear  problems  with  prescribed  boundaries.  Finally, 
because  the  wetted  part  of  wetland  may  change  significantly  with  variations  of  wetland  water 
stage,  the  finite-element  mesh  must  be  adjusted  horizontally  and  vertically  between  time 
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steps. 

Before  the  pumping  test  began,  the  SV5  wetland  water  surface  was  about  60  meters 
in  diameter.  However,  after  pumping  stopped  it  was  less  than  30  meters.  Thus,  the  horizontal 
inundated  area  changed  from  2750  m2  to  750  m2.  To  model  the  influence  of  a  changing 
horizontal  free  surface  required  that  the  finite -element  mesh  be  adjusted  horizontally.  The 
zone  of  mesh  adjustment  would  be  within  the  peat  layer  and  include  an  upper  boundary 
representing  the  free  surface  or  the  surface  water  stages  and  a  lower  boundary  corresponding 
to  the  vertical  limit  of  the  peat  layer.  Adjusting  the  mesh  horizontally  keeps  the  same  number 
of  nodes  on  the  boundary.  Thus,  the  physical  domain  of  the  mesh  may  change  horizontally 
and  vertically.  However,  the  calculation  domain  does  not.  After  each  simulation,  information 
defining  the  location  of  every  node  is  updated. 

Since  the  wetland  boundary  was  assumed  to  vary  over  time,  the  total  head  was  taken 
as  a  boundary  condition.  Thus,  at  any  time,  the  wetland  boundary  can  be  viewed  as  an 
equipotential  line.  Water  moves  vertically  to  and  from  the  equipotential  line  represented  by 
the  free  surface  boundary.  The  vertical  hydraulic  gradient  at  the  wetland  boundary  times  the 
product  of  the  wetland  area  and  hydraulic  conductivity  equals  the  volume  of  water  exchanged 
between  the  wetland  and  the  aquifer. 
3.6.6  Model  Calibration  and  Results 

Because  the  aquifer  underlying  SV5  is  composed  of  sand,  a  specific  storativity  of 
S ,  =  10"5  m"1  was  used  (Bear  1979).  The  saturated  water  content  and  residual  water  content 
in  the  peat  layer  were  0S  =  0.7  and  0r  =  0.6.  In  the  aquifer,  0S  =  0.34  and  0r  =  0.18,  based  on 
Walser  (1998).  Water  stage  level  observation  data  and  groundwater  monitoring  data  collect 
during  the  pump  test  and  many  hours  after  pumping  stopped  were  used  to  calibrate  the 
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numerical  model. 

In  Krabbenhoft's  (1996)  interaction  model  between  groundwater  and  lakes,  an 
anisotropy  ratio  of  100  was  obtained  after  several  trial  and  error  simulations.  Also  an 
anisotropy  ratio  of  horizontal  to  vertical  hydraulic  conductivity  of  10  was  obtained  by  Winter 
(1978)  in  his  numerical  simulation  of  steady  state  three-dimensional  groundwater  flow  model 
near  lakes. 

In  this  finite-element  model,  it  was  assumed  that  the  hydraulic  conductivity  of  the 
peat  layer  was  different  from  the  hydraulic  conductivity  of  the  underlying  aquifer,  and  the 
aquifer  was  also  assumed  anisotropic.  Furthermore,  it  was  assumed  that  the  hydraulic 
conductivities  and  anisotropies  of  both  layers  were  spatially  uniform  within  each  layer 
throughout  the  aquifer  zone.  By  adjusting  the  anisotropy  ratio  (Kh/Kv)  and  checking  for  the 
best  agreement  between  modeled  and  observed  data,  it  was  determined  that  the  best  results 
were  obtained  when  Kh/Kv=  1.0,  Kxx=  K^  =  8.7  m/day,  and  K^,  =  8.7  m/day  in  the  aquifer, 
while  Kp  =  0.2  m/day  in  the  peat  layer.  This  number  is  different  from  the  value  of  0.03m/day 
obtained  by  Switt  (1998)  and  0.087m/day  obtained  by  Wise  (2000).  The  possible  difference 
may  be  caused  by  several  reasons.  First,  the  peat  layer  depth  used  in  the  numerical  model  is 
not  a  constant.  When  the  water  stage  in  the  wetland  is  lowered,  the  peat  layer  depth  in  the 
inundated  area  is  much  greater  than  the  average  value  of  0.25  m  used  in  the  analytical  model. 
The  variation  of  groundwater  head  immediately  beneath  the  peat  layer  may  not  be  the  same 
value  estimated  using  a  groundwater  head  variation  at  a  depth  of  6  meters  beneath  of  the 
wetland.  Thus,  the  second  reason  may  be  that  the  head  distribution  beneath  the  wetland 
obtained  from  the  numerical  model  is  different  from  the  values  used  by  the  analytical  model. 

Since  the  aquifer  beneath  the  SV5  wetland  is  composed  of  sand,  the  specific 
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storativity  is  extremely  small.  The  heads  along  the  depth  are  very  different.  Figure  3-27 
shows  the  head  distribution  beneath  the  wetland.  Figure  3-28  shows  the  head  variation  with 
time  at  different  depths  beneath  the  wetland.  The  piezometric  head  changes  relatively  quickly 
in  the  deep  zone,  but  the  head  does  not  change  much  near  the  wetland.  Figure  3.29  shows  the 
variation  of  the  inundated  area  of  the  wetland.  The  area  changes  from  2,600  m2  to  700  m2just 
after  pumping.  Figure  3-30  shows  the  variation  of  the  interaction  discharge.  The  discharge 
supplied  from  the  aquifer  to  the  wetland  increases  quickly  and  then  decreases  slowly  after 
pumping  was  stopped.  Figure  3-31  shows  excellent  agreement  between  the  simulated  results 
given  by  the  saturated  flow  model  and  observed  groundwater  levels.  Table  3-4  to  Table  3-8 
show  the  modeled  data  used  to  calculate  the  interaction  discharges  at  5  different  times. 

Table  3-8  shows  the  sensitivities  of  the  head  relative  to  changes  in  the  hydraulic 
conductivities  in  the  pear  layer  and  the  aquifer  and  to  changes  in  the  specific  storativity.  The 
hydraulic  conductivities  in  both  the  peat  layer  and  the  aquifer  play  the  most  important  part 
in  the  change  of  head  distribution  and  interaction  discharge.  The  sensitivity  of  the  specific 
storativity  is  relatively  small. 
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Figure  3-27.  Vertical  cross-section  view  of  three-dimensional  finite-element 
mesh  and  head  from  the  saturated  flow  model  at  a  time  of  87  hours. 
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Figure  3-28.  Piezometric  head  at  positions  beneath  the  wetland  with  time 
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Figure  3-29.  Variation  in  the  inundated  area  of  wetland  during  simulation 
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Figure  3-30.  Discharge  between  wetland  and  aquifer  with  respect  to  time 
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Figure  3-31.  Comparison  between  calculated  head  and  observed  head 


Table  3-4.  Data  for  calculating  the  discharge  to  the  wetland  from  the  aquifer  at  10  hours 


Elements 
of  Peat  Layer 

Area 

(m2) 

Head  in  Peat 
Layer      (m) 

Head  on  Peat 
Surface       (m) 

Distance  of  Two 
Nodes            (m) 

1 

6.95828 

15.7810 

15.7223 

0.301929 

2 

20.8748 

15.7813 

15.7223 

0.305786 

3 

34.7914 

15.7819 

15.7223 

0.306822 

4 

48.7075 

15.7827 

15.7223 

0.303675 

5 

62.6240 

15.7843 

15.7223 

0.297326 

6 

76.5414 

15.7868 

15.7223 

0.289907 

7 

90.4590 

15.7898 

15.7223 

0.281774 

8 

104.377 

15.7936 

15.7223 

0.271725 

9 

118.296 

15.7979 

15.7223 

0.260863 

10 

132.215 

15.8028 

15.7223 

0.248279 

Q  to  wetland  at  10  hr 

1.48466  mVhr 
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Table  3-5.  Data  for  calculating  the  discharge  to  the  wetland  from  the  aquifer  at  20  hours 


Elements 
of  Peat  Layer 

Area 
(m2) 

Head  in  Peat 

Layer    (m) 

Head  on  Peat 
Surface     (m) 

Distance  of  Two 
Nodes         (m) 

1 

8.36449 

15.7953 

15.7427 

0.302115 

2 

25.0935 

15.7956 

15.7427 

0.306344 

3 

41.8223 

15.7960 

15.7427 

0.306027 

4 

58.5508 

15.7964 

15.7427 

0.302183 

5 

75.2802 

15.7982 

15.7427 

0.294043 

6 

92.0102 

15.8003 

15.7427 

0.287151 

7 

108.741 

15.8036 

15.7427 

0.275314 

8 

125.472 

15.8070 

15.7427 

0.264839 

9 

142.204 

15.8114 

15.7427 

0.250437 

10 

158.937 

15.8162 

15.7427 

0.236900 

Q  to  wetland  at  20hr 

1.65871  m3/hr 

Table  3-6.  Data  for  calculating  the  discharge  to  the  wetland  from  the  aquifer  at  40  hours 


Elements 
of  Peat  Layer 

Area 
(m2) 

Head  in  Peat 
Layer     (m) 

Head  on  Peat 
Surface    (m) 

Distance  of  Two 
Nodes       (m) 

1 

10.9288 

15.8222 

15.7797 

0.302417 

2 

32.7864 

15.8225 

15.7797 

0.307252 

3 

54.6436 

15.8225 

15.7797 

0.304707 

4 

76.5010 

15.8229 

15.7797 

0.298378 

5 

98.3597 

15.8243 

15.7797 

0.289097 

6 

120.219 

15.8264 

15.7797 

0.277635 

7 

142.080 

15.8288 

15.7797 

0.266032 

8 

163.941 

15.8321 

15.7797 

0.249581 

9 

185.804 

15.8359 

15.7797 

0.233460 

10 

207.669 

15.8406 

15.7797 

0.216548 

Q  to  wetland  at  40hr 

1.86649  m3/hr 
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Table  3-7.  Data  for  calculating  the  discharge  to  the  wetland  from  the  aquifer  at  60  hours 


Elements 
of  Peat  Layer 

Area 
(m2) 

Head  in  Peat 

Layer    (m) 

Head  on  Peat 
Surface    (m) 

Distance  of  Two 
Nodes      (m) 

1 

13.1219 

15.8463 

15.8115 

0.302649 

2 

39.3657 

15.8466 

15.8115 

0.307800 

3 

65.6089 

15.8464 

15.8115 

0.303651 

4 

91.8530 

15.8469 

15.8115 

0.295196 

5 

118.098 

15.8480 

15.8115 

0.285844 

6 

144.345 

15.8498 

15.8115 

0.271268 

7 

170.593 

15.8519 

15.8115 

0.255422 

8 

196.843 

15.8544 

15.8115 

0.239015 

9 

223.095 

15.8579 

15.8115 

0.219572 

10 

249.349 

15.8626 

15.8115 

0.199788 

Q  to  wetland  at  60hr 

1.94458  rnVhr 

Table  3-8.  Data  for  calculating  the  discharge  to  the  wetland  from  the  aquifer  at  80  hours 


Elements 
of  Peat  Layer 

Area 

(m2) 

Head  in  Peat 
Layer    (m) 

Head  on  Peat 
Surface    (m) 

Distance  of  Two 
Nodes      (m) 

1 

15.0334 

15.8681 

15.8391 

0.302835 

2 

45.1003 

15.8683 

15.8391 

0.307522 

3 

75.1664 

15.8680 

15.8391 

0.302700 

4 

105.234 

15.8685 

15.8391 

0.292596 

5 

135.303 

15.8695 

15.8391 

0.280887 

6 

165.374 

15.8709 

15.8391 

0.266576 

7 

195.446 

15.8727 

15.8391 

0.247805 

8 

225.521 

15.8749 

15.8391 

0.228349 

9 

255.599 

15.8777 

15.8391 

0.208494 

10 

285.679 

15.8819 

15.8391 

0.186883 

Q  to  wetland  at  80hr 

1.93622  m'/hr 
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Table  3-9.  Sensitivities  calculated  by  saturated  groundwater  flow  model 


Parameter  (PR) 

Percentage  varied 

APR 

llh,-h2ll  (m) 

llh,-h2ll/APR 

Ka(m/day) 

-0.10 

-0.85 

0.9924E-02 

-1.1675E-02 

Ka(m/day) 

+0.10 

+0.85 

0.9055E-02 

1.0653E-02 

Kp(m/day) 

-0.10 

-0.02 

0.6410E-02 

-0.3205 

Kp(m/day) 

+0.10 

+0.02 

0.5564E-02 

0.2782 

S0 

-0.50 

-  0.000005 

0.2068E-04 

-4.136 

S0 

+0.50 

+0.000005 

0.1265E-04 

2.53 

3.7  Comparison  with  Results  of  Wise's  (2000)  and  others 

Parameter  values  obtained  by  Wise  (2000),  1^=  0.087  (m/day),  by  Walser  (1998), 
kSand  =  0-75  (m/day),  kpeat  =  0.03  (m/day)  for  SV5  were  put  in  the  three-dimension  finite- 
element  model  in  this  dissertation,  and  the  simulated  results  were  plotted  on  the  Figure  3-32. 
When  the  specific  storage  was  set  on  0.01  (1/m),  both  the  simulated  wetland  water  table  and 
groundwater  are  pretty  close  to  the  observed  conditions.  The  simulated  wetland  water  table 
and  the  groundwater  table  corresponding  to  the  specific  storage  of  0.001  (1/m)  are  much 
different  from  the  observed  results. 

There  are  two  main  differences  between  Walser' s  numerical  simulation  and  the 
models  in  this  dissertation.  First,  unlike  a  coupling  boundary  condition  used  in  this 
dissertation,  the  fixed  observed  water  table  process  over  time  was  used  in  Walser' s  numerical 
simulation.  The  other  difference  is  that  a  fixed  three-dimension  mesh  was  used  in  Walser' s 
numerical  simulation. 

When  kaqu  at  8m/day  and  k^at  0.2m/day  were  used ,  and  specific  storage  was  varied, 
the  simulated  results  are  plotted  on  the  Figure  3-33.  The  simulated  wetland  water  tables  for 
both  specific  storage  values  are  very  close  to  the  observed  wetland  water  table  results.  The 
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groundwater  table  for  the  specific  storage  of  0.01(l/m),  however,  is  above  the  observed 
results. 
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Figure  3-32.  Comparison  results  using  Walser's  parameter  values. 
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Figure  3-33.  3-D  numerical  model  simulated  results  over  specific  storage 
Table  3-10  shows  the  hydraulic  conductivities  in  sand  layer  and  peat  layer.  The 
variations  are  rather  large. 
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Table  3-10.  Hydraulic  Conductivity  values  for  sand  and  peat  layers 


k  in  sand  layer 
(m/day) 

k  in  peat  layer 
(m/day) 

Swittetal.(1998) 

7-12 

0.03 

Slug  test 

2 

Rycroft  et  al.  (1975) 

10"7-102 

Wise  et  al.  (000) 

0.083 

Walser  Numerical  Model  (1998) 

0.75 

0.03/0.3 

Models  in  this  dissertation 

8 

0.2 

3.8  Summary 

A  three-dimensional  finite -element  numerical  saturated  groundwater  flow  model 
using  deformable  mesh  and  real-time  wetland-aquifer  coupling  boundary  is  developed  in  this 
chapter.  Considering  the  area  of  the  water  surface  changes  quickly  because  most  wetlands 
have  a  shallow  bank  slope,  the  model  has  the  capability  of  tracking  the  interaction  the 
interaction  boundary  between  groundwater  and  a  wetland.  The  model  is  verified  using  two 
analytical  solutions,  one  is  Theis  solution,  the  other  is  the  analytical  solution  found  in 
Chapter  2  of  this  dissertation,  which  describes  groundwater  flow  around  a  circular  wetland. 
The  model  is  further  used  to  simulate  the  interaction  problem  of  groundwater  and  the  SV5 
wetland.  The  pumping  from  the  SV5  wetland  and  the  recovery  of  the  SV5  wetland  surface 
afterwards  are  closely  simulated  by  this  model,  the  interaction  boundaries  are  effectively 
traced  by  the  deformable  mesh. 


CHAPTER  4 

INVERSE  MODEL  OF  THREE-DIMENSIONAL  FINITE-ELEMENT  METHOD 

SATURATED  GROUNDWATER  MODEL  IN  SEARCHING 

FOR  THE  HYDRAULIC  CONDUCTIVITY  OF  THE  PEAT  LAYER 


4.1  Introduction 
4.1.1  Inverse  Problem 

During  the  past  twenty  years,  mathematical  modeling  has  been  extensively  used  in 
the  study  of  groundwater  problems.  However,  mathematical  models  require  a  number  of 
parameters  to  characterize  the  aquifer.  Usually,  prior  knowledge  of  the  values  of  these 
parameters  is  limited,  and  they  must  be  quantified  by  calibration  which  requires  observations 
of  hydraulic  heads  or  solute  concentrations. 

Given  a  parametric  model  of  the  physical  system  and  values  of  the  model  parameters, 
the  prediction  of  system  output  (response)  to  any  input  is  referred  to  as  the  direct  problem, 
or  the  forward  problem  (simulation).  The  forward  problem  predicts  unknown  system  states 
by  solving  appropriate  governing  equations. 

Conversely,  the  estimation  of  model  parameters  given  the  parameteric  system  model 
and  an  input-output  relationship  is  known  as  the  inverse  problem  (calibration).  The  typical 
inverse  problem  involves  estimating  the  values  of  model  parameters,  solving  the  direct 
problem  with  this  set  of  parameters  to  calculate  the  predicted  system  response,  and  then 
adjusting  model  parameters  until  deviations  between  observed  and  model  predicted 
responses  are  in  some  specified  sense  minimized  (Sun  1994). 
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Trial  and  error  and  graphical  matching  techniques  are  the  most  basic  and  traditional 
methods  to  solve  inverse  problem.  They  are  widely  used  because  they  are  simple  to 
implement.  Because  the  number  and  combinations  of  parameter  adjustments  are  not 
bounded,  "trial  and  error"  is  rather  flexible  but  also  time-consuming.  Additionally,  the 
solution  is  strongly  dependent  on  the  skill  of  the  practitioner  (Keidser  1991). 
4.1.2  Parameterization 

Most  groundwater  inverse  algorithms  adopt  either  a  blocked  or  a  geostatistical 
(random  field)  description  of  spatial  variability.  The  first  approach  divides  the  region  of 
interest  into  a  number  of  discrete  blocks  or  zones  that  are  believed  to  correspond  to  distinct 
geological  units  (Yeh  1986,  Carrera  and  Neuman  1986a,b,c,  and  Cooley 
1977,1979,1982,1983).  Each  block  is  characterized  by  a  set  of  spatially  uniform 
hydrogeologic  properties  which  are  treated  as  parameters  in  an  appropriate  inverse  problem. 
Thus,  the  unknown  distributed  paramters  can  be  represented  by  a  number  of  constants.  The 
weakness  of  this  method  is  obvious,  i.e.,  although  we  have  some  of  the  geological 
information  is  available,  it  is  still  difficult  to  divide  the  flow  region  into  a  limited  number  of 
zones.  If  the  zonation  pattern  is  incorrect,  the  identified  paramter  values  also  will  be  incorrect 
(Sun  and  Yeh  1985). 

The  geostatistical  alternative  views  the  properties  of  interest  as  stationary  random 
fields  that  vary  relatively  smoothly  over  space  (Hoeksma  and  Kitanidis  1984,  Dagan 
1985,Yeh  1996,  and  Zhang  1997,  Michalak  2004).  In  the  geostatistical  interpolation  method 
of  kriging,  the  unknown  parameter  is  regarded  as  a  stochastic  field  described  by  some 
statistical  parameters,  such  as  the  mean,  variance,  and  correlation  distance.  Then  the 
identified  parameter  vector  consists  of  only  a  few  statistical  parameters  if  the  aquifer  is 
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basically  homogeneous.  The  problem  of  overparameterization  may  be  avoided,  and  the 
inverse  solution  obtained  by  the  maximum  likelihood  estimate  and  cokriging  is  generally 
stable.  The  correctness  of  the  inverse  solution,  however,  depends  on  the  correctness  of  some 
statistical  assumptions  and  the  structure  of  covariance  functions. 

Although  the  two  approaches  are  based  on  different  parameterizations,  they  both  treat 
hydrogeologic  properties  as  spatial  functions.  This  suggests  that  it  should  be  possible  to 
formulate  a  general  inverse  theory  which  encompasses  both  the  blocked  and  geostatistical 
alternatives,  as  well  as  hybird  methods  which  lie  between  these  extremes. 
4.1.3  Parameter  Identification  Methods 

Many  methods  are  available  for  the  solution  of  the  common  groundwater  modeling 
inverse  problem  whereby  the  conductivity  or  transmissivity  of  a  porous  medium  are 
estimated  using  available  observations  of  head  or  pressure  and  other  information  (Distefano 
1975,  Chang  1976,  Loaiciga  1987,  Mishra  1989,  Knopman  1989,  Keidser  1991,  Sun  1995, 
Yeh  1996,  Zhang  1997,  Lamorey  1998,  Anderman  1999,  Michalak  2004).  Neuman  (1973) 
classified  inverse  modeling  solution  techniques  as  either  "direct"  or  "indirect". 

The  direct  method  or  equation  error  method  is  used  when  head  variations  and 
derivatives  (usually  estimated)  are  known  over  the  entire  flow  region.  If  the  measurement 
and  model  errors  are  negligible,  the  original  governing  equation  becomes  a  linear  first-order 
partial  differential  equation  of  the  hyperbolic  type  in  terms  of  the  unknown  parameters.  With 
the  aid  of  boundary  conditions  and  flow  data,  a  direct  solution  for  the  unknown  parameters 
may  be  possible.  In  practice,  observation  wells  are  sparsely  distributed  in  the  flow  region  and 
only  a  limited  number  of  observation  wells  are  available.  To  formulate  the  inverse  problem 
by  the  equation  error  criterion,  missing  data  have  to  be  estimated  by  interpolation.  The 
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interpolation  data  contain  errors  resulting  from  the  interpolation,  and,  thus,  cause  additional 
errors  in  the  results  of  parameter  identification. 

The  solution  approach  of  the  indirect  method  to  solve  inverse  problems  is  to 
formulate  the  inverse  problem  in  term  of  an  optimization  framework,  based  on  the  most 
common  output  error  criterion,  i.e.,  output  error  least  squares  criterion.  Using  this  approach, 
a  set  of  parameters  is  identified  which  best  describes  observed  system  states. 

With  inverse  modeling,  the  structure  and  parameters  of  a  model  are  adjusted 
simultaneously  or  sequentially  so  as  to  obtain  a  parameter  input  system  and  output  relations 
that  characterize  observed  excitation-response  functions  of  the  real  system.  Since  the  model 
structure  for  the  wetland/groundwater  system  is  assumed  to  follow  Eq.5-1-5-3,  the  problem 
of  determining  the  hydraulic  conductivity  of  the  peat  layer  from  the  observed  water  stage  and 
total  head  in  the  aquifer  may  be  taken  as  a  typical  inverse  problem.  The  investigation  of  the 
inverse  problem  based  on  the  three-dimensional  saturated  groundwater  flow  model  and  the 
adjoint  state  method  is  presented  in  this  chapter. 
4.1.4  Computation  of  Sensitivity  Coefficients 

Sensitivity  coefficients,  i.e.,  the  partial  derivatives  of  head  with  respect  to  each  of  the 
parameters,  play  an  important  role  in  the  solution  of  the  inverse  problem.  In  the  Gauss- 
Newton  algorithm,  elements  of  the  Jacobian  matrix  are  represented  by  the  sensitivity 
coefficients,  dh/dT„  i=l,....,M,  1=1,.. .,L.  If  h  is  the  head  vector,  the  sensitivity  coefficients 
are  dh/dT„  1=1,. ..,L.  There  are  three  methods  used  in  calculation  of  sensitivity  coefficients. 
There  are  the  influence  coefficient  method,  sensitivity  equation  method,  and  variational 
method.  The  influence  coefficient  method  determines  the  sensitivity  coefficients  by 
perturbing  each  parameter  once  at  a  time. 
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In  the  sensitivity  equation  method,  a  set  of  sensitivity  equations  is  obtained  by  taking 
the  partial  derivatives  with  respect  to  each  parameter  in  the  governing  equation  and  in  the 
initial  and  boundary  conditions.  The  sensitivity  coefficients  are  obtained  by  solving  the  new 
governing  equations. 

The  variational  method  was  first  used  for  solving  the  inverse  problem  of  parameter 
identification  by  Jacquard  and  Jain  (1965)  and  then  by  Carter  et  al.  (1974,  1982)  associated 
with  finite-difference  schemes.  Sun  and  Yeh  (1985)  extended  the  method  to  the  case  of  a 
finite-element  scheme.  The  sensitivity  coefficients  can  be  computed  using  the  following 
equation: 


dhm        '  V 


3lf 

0   Qt 


dq{xy,t-x) 
{  dt 


•VACW)  cbufydx  (4"1) 


j=l,2,...,N0  i=l,2,...,Nn 

where  Q  j  is  the  exclusive  subdomain  of  node  i;  V  is  the  gradient  operator;  h(x,y,t)  is  the 

solution  of  the  governing  equation;  N0  is  the  number  of  observation  wells;  Nn  is  the  total 

number  of  nodes  used  in  the  numerical  solution;  and  q(x,y,t)  is  the  time  derivative  of  q(x,y,t), 

which  is  the  solution  of  the  set  of  adjoint  equations. 

4.2  Adjoint  Problem  for  Three-Dimensional  Finite-Element 
Saturated  Groundwater  Model 

The  focus  of  this  section  is  the  formulation  of  the  adjoint  problem  which  defines 

sensitivities  to  be  used  in  determining  the  hydraulic  characteristics  of  the  peat  layer  and  the 

underlying  aquifer.  It  is  assumed  that  the  thickness  of  the  peat  and  aquifer  are  known.  As 
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discussed  in  Chapter  3,  three  dimensional  groundwater  flow  in  an  unconfined  aquifer  is 
governed  by: 


dt  J 


(4-2) 


subject  to  the  initial  and  boundary  conditions: 


h(x,y,z,0)  =f0  (xyj)€Q 


(4-3) 


*lr,  =/i» 


KVh-nL  =/,  , 


(4-4) 


f0  is  the  initial  head;  f,  is  the  specific  head  condition;  f,  is  the  specified  flux  condition;  and 
tf  is  the  time  duration.  The  boundary  of  volume  Q  shown  in  Figure  4-1  is  denoted  by  T 
including  wetland  head  boundary  (I\),  no  flow  boundary  (IT1,),  and  lateral  head  boundary 

(r2). 


A 


Z  _ 


£~ ft    ^ 


Peat  Layer- 

Groundwater  Aquifer 


HI 


AB,  EG,  FH  or  CD  are  head  boundaries  T2 
GH  or  EC,  DF  are  no  flow  boundaries  T  , 


Figure  4-1.  Sketch  showing  setup  of  the  saturated  groundwater 
inverse  problem  around  a  wetland 
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The  hydraulic  conductivities  of  the  peat  layer  and  the  underlying  aquifer  are  the 
unknown  parameters.  Any  variation  SK  of  K  where  K=hydraulic  conductivity  must  cause  a 
variation  6h  of  head  h,  because  h  is  dependent  on  K  through  the  governing  equation  and 
subsidiary  conditions.  Taking  the  variation  of  Eq.  4-2  and  subsidiary  conditions,  the 
following  is  obtained: 


V-(5KVh)+V{KV6h)=Ss^Q         (x,y,z)eQ,       0<t<tf  (4.5) 


subject  to  the  conditions: 

6h(x,y,z)\t_0  =  0  OwOeQ  (4.6) 

8/2|p2  =  0,  [6KVh+KV6h]-n\T   =0  ,  (4.7) 

Eq.  4-5-Eq.  4-7  represent  the  variational  problem  of  the  primary  problem  in  Eq.  4-2.  The 
variational  problem  relates  the  variation  6k  and  oh.  Assuming  h  is  an  arbitrary  function 
having  continuous  second-order  space  derivatives  on  Q  and  first-order  derivatives  in  [0,  tf], 
then  multiplying  Eq.  4-5  by  T  and  integrating  the  result  over  time-space  domains  yields: 


(SsV6h)\'{dQ- 

0    Q 


Ss6h—  +  WVidkVh)  +  ¥V-0tV6/2) 


dQdt  -  0  (4-8a) 
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If  we  define: 


V(x,y,z,tf)  =  0  (x^)eQ 


(4-8b) 


then  the  first  term  of  Eq.  4-8a  vanishes  because  of  condition  4-7.  Thus,  Eq.  4-8b  reduces  to: 


Ss5h—  +  YV-(6£V/z)  +  VVikVbh) 


dQdt  =  0 


o  a 


(4-9) 


where: 


/ 


YV-(5A:V/0dQ  -  |  y(6kVhynds  -  |  6kVhWdQ 

a  s 


(4-10) 


and: 


¥VikVbh)dQ  =  I  ViWbhynds  -  \  Wdh-VYdQ 
a 


(4-11) 


-     Wbh-WdQ  =  J  6h\yikVW)]dQ  -   J  bhkVU-nds 


(4-12) 


Using  identities  of  Eq.  4-10,  Eq.  4-1 1  and  Eq.  4-12,  Eq.  4-9  can  be  written  as: 
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aw 

Sfih—  +  6hV7'(WT)]  -fiAV/rVP 

*      dt 


dQ 


dt 


ybkVh-nds+    ^kVbh-nds-     bhWS-nds 


dt=0 


(4-13) 


Using  conditions  4-7  and  defining: 


¥1      =0 


kVY-n\r  =  0 


then  the  second  term  of  Eq.  4-13  is  equal  to  zero.  Eq.  4-13  produces: 


s  dt         K        ' 


dhdQdt  -        {V¥-Vh)dkdQdt  =  0 


o  o 


o  o. 


(4-14) 


In  order  to  identify  the  hydraulic  conductivity  k(x,y,z)  that  would  lead  to  the  best  prediction 
of  observed  heads,  a  performance  criterion  must  be  defined.  A  general  form  of  the 
performance  function  may  be  written  as: 


E(h,k)  =        j{h,k;x,y,z,t)dQdt  =  0 
o  a 


(4-15) 
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Taking  the  variation  of  E(h,  k)  yields: 


dE(h,k)  =  I  I  [  ^bk+^bh\dQdt  = 
]]{dk       dh     ) 

o  a 


(4-16) 


The  summation  of  Eq.  4-16  and  Eq.  4-14  gives: 


6E  = 


dt     s  dt 


(AVP) 


0    Q 


o  a 


\ 


dkdQdt  (4"17) 


If  T(x,y,z,t)  is  selected  to  satisfy  the  following  PDE: 


S  —  +V-(A:V,P)--^  =  0 

*  ar  '   a/2 


(4-18) 


subject  to  conditions: 


Y(x,>>,z,f,)  =  o 


WL  =  0,  *VY-#iL  =0 


then,  the  first  integral  on  the  right  hand  side  of  Eq.  4-17  related  to  unknown  8h  will  vanish. 
Thus: 


6E 


(v 

dk 


+  VT-V/Z 
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dk  dQ  dt 


o  a 


From  Eq.  4-19,  the  following  partial  derivative  can  be  defined: 


5E 

dk 


M. 

[dk 


+  VY-V/? 


dQdt 


o  o. 


By  introducing  the  following  time  transformation: 


(4-19) 


(4-20) 


T-tf-t 


(4-21) 


the  adjoint  problem  can  be  expressed  as 


5  — -V-(*VP)  +  -^  =  0 

*  dx  dh 


(4-22) 


subject  to: 


*l,=o=°.  (*.>oea 


xli\=o   u 


*VY-ii|r  =0 


(4-22a) 


where  Y  =  Y(x,y,z,tf  -t).  Solving  each  equation  once,  all  components  of  the  gradient  VE(k) 
can  be  obtained. 
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4.3  Searching  Objective  Function 

To  find  the  best  parameter  value,  the  model  structure  and  model  parameters  are 
solved  simultaneously  or  sequentially  so  as  to  make  the  input-output  relation  of  the  model 
fit  any  observed  excitation-response  relationship  of  the  real  system.  Since  the  model  structure 
of  the  problem  in  this  study  is  determined,  the  problem  of  determining  the  hydraulic 
conductivity  of  the  peat  layer  and  the  underlying  aquifer  from  the  observed  water  stage  and 
total  head  in  the  aquifer  may  be  taken  as  a  typical  inverse  problem: 


E(Kest)=minE(K) 


(4-23) 


E(k)-EKrKJ  (4-24) 


where: 


L 

I 

The  traditional  method  to  solve  this  kind  of  problems  is  using  the  trial  and  error  method 
based  on  the  output  least  squares  criterion  (Eq.  4-23). 

4.4  Adjoint  State  Controlling  Equations 
In  this  study  the  adjoint  state  method  was  used  to  solve  the  inverse  problem.  From 
the  governing  equation  (Eq.  4-2)  and  the  associated  boundary  conditions,  the  adjoint 
formulation  of  the  primary  problem  may  be  found  using  Green's  first  and  second  theorems 
and  other  transforms: 


Ss  —  -V-(fcVY)  +  -^  =  0 

dx  dh  (4.25) 
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subject  to: 


Yl    n=0,    (x,y)eQ 
It-o     »    v  ^  (426) 


and 


Tlryo=°  *VT-»|r=0  (4.27) 


where  h  is  the  total  head;  T  =  Y  (x.y^-T)  is  the  adjoint  state  of  h;  tf  is  the  time  duration  of 
the  simulation  criterion  and  f(h,k;x,y,t)  expresses  the  'goodness'  of  fit  between  the  predicted 
head,  h,  and  the  observed  head,  hobs  expressed  through  the  ordinary  least  squares  (OLS): 


I 

i  1 


M*VA0  =  E  I  h(k)-hob\\6(x-xobs)6(y-yobs)6(z-zobs)6(t-tp  (4-28) 


where  (  xobs,  yobs,  zobs)  is  the  location  of  the  observation  well;  and  tj  (j  =1,2, "-J)  are  the 
observation  times.  Integrating,  Eq.  4-28  over  time  space  domains  produces: 


E(h,k)=j     j     Y,f{h,k-x,y,z,t)dQdt  (4-29) 

Jo    «/(Q)  j=l 


By  solving  the  governing  Equation  4-2  and  the  adjoint  Equation  4-25  once,  we  may  obtain 
the  following  functional  partial  derivative  required  to  solve  Equation  4-23  is  obatined: 
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bE 
bk 

'0    «/(Q) 


dk 


dQdt 


(4-30) 


4.5  Searching  Method  and  Stopping  Criterion 

Since  the  gradient  VE(k)  is  known,  many  searching  methods  including  the  gradient 
method,  quasi-Newton  method,  and  the  Fletcher-Reeves  method  could  be  used  to  solve  the 
inverse  problem  4-23  (Sun  1985).  Two  stopping  criteria  used  here  are: 


■w.i<»,  (4.31) 


and: 


\\E(kn)-E(kn^\\<e2  (4-32) 

Using  the  adjoint  method  and  the  Fletcher-Reeves  algorithm  to  update  the  search  sequence, 
a  highly  efficient  numerical  scheme  for  solving  the  inverse  problem  can  be  formulated. 

4.6  Field  Application  on  SV5 
The  inverse  model  was  used  to  estimate  the  hydraulic  conductivity  of  the  peat  layer 
underlying  the  SV5  wetland  (Wise  et  al.  2000).  Since  the  piezometric  head  at  a  depth  of  5 
meters  beneath  the  wetland  was  observed,  f  is  defined  as: 

AMW*0-£  |  Hk)  -hoblb(x-xobs)b(y-yob)b(z-zob)b{t-t)  (4-33) 

7=1 

Water-level  data  were  measured  every  hour  for  a  total  duration,  tf,  of  87  hours.  An  example 
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of  how  the  hydraulic  conductivity,  K,  is  updated  between  iterations  is  as  follows.  Assuming 
an  initial  estimate  of  the  vertical  hydraulic  conductivity  K  of  the  peat  layer  was 
K,=1.00m/day,  then  8E/aK  =  -0.0125  by  Eq.4-29  and  E(K,)=  0.48998.  Because  it  is  not 
possible  to  match  the  modeling  result  with  the  observed  result  exactly,  the  final  E(k)  is  not 
zero.  The  sensitivity  was  used  to  decide  the  direction  of  searching.  Assuming  E(K)  is 
supposed  to  be  the  smallest  value,  the  updated  value  of  K=K2  is  assumed  to  be  0.6  m/day  as 
generated  from  the  formula: 


E(K2)=E{Kx)  +  ^\     {K.-K,)  (4-34) 

oK      1 


Based  on  the  new  K  =  K2  =  0.6  m/day,  the  updated  gradient  dE/dK  =  0.0188,  and  E(K2)= 
0.41467.  Then  a  new  K  is  assumed  until  E(K)  is  minimized.  The  final  K  is  0.2  m/day. 
Figure  4-2  and  Table  4-1  are  the  calculated  values  from  the  adjoint  method  model.  In  this 
model,  the  observed  groundwater  data  were  used  in  the  output  least  squares  criterion.  The 
observed  water  stage  in  the  wetland  was  not  used  in  the  model.  This  is  the  major  reason  that 
the  final  result  from  the  inverse  model  is  a  little  bit  different  from  the  result  obtained  in 
Chapter  3.  The  saturated  flow  model  in  Chapter  3  uses  both  the  observation  of  groundwater 
head  beneath  the  wetland  and  the  surface  water  stage  of  wetland  to  calibrate  the  model.  In 
the  inverse  model,  however,  it  is  hard  to  use  the  surface-water  observation  result,  since  it  is 
a  boundary  of  the  model. 

Because  the  wetted  area  of  surface  water  associated  with  the  wetland  changes  with 
time,  the  space  location  of  each  finite-element  point  tied  to  the  location  of  the  water  surface 
must  be  saved  before  the  adjoint  state  equations  are  solved.  Thus,  for  each  assumed  hydraulic 
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conductivity  in  the  peat  layer,  the  forward  problem,  and  adjoint  state  problem  must  be 
solved  first,  and  then  <3E/dk  for  this  assumed  value  of  K  could  be  found.  Figure  4-3  shows 
the  procedure  for  solving  an  inverse  problem  using  adjoint  method. 
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Figure  4-2.  Relation  between  error  and  peat  hydraulic  conductivity 


Table  4-1.  Error  and  sensitivity  based  on  the  adjoint  method  model 


K  (m/day) 

E(k) 

6E/aK 

1 

1.0 

10.5698 

0.26965 

2 

0.6 

8.94509 

0.40555 

3 

0.4 

7.17941 

1.03330 

4 

0.2 

3.47742 

0.31495 

5 

0.1 

4.82392 

3.72980 

6 

0.08 

5.59932 

3.14520 
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y. 

Model  groundwater  flow 
based  on  the  assumed  parameters  K„ 


^ 


Solve  the  adjoint  state  governing  equations 


±. 


Calculate  dE/dK  based  on  the  distributions 
of  groundwater  head  and  adjoint  state 


Calculate  E  and  new  K^, 


No 


Yes 


V_ 


END 


A 


Figure  4-3.  Flowchart  illustrating  the  iteration  procedure  of  inverse  problem 

4.8  Summary 
The  inverse  problem  for  three-dimensional  transient  groundwater  flow  in  a  phreatic 
aquifer  is  more  complicated,  because  the  free  surface  changes  with  time  such  that  the 
simulation  zone  is  not  stationary.  In  order  to  solve  the  adjoint  state  governing  equations,  the 
position  of  each  node  must  be  saved  at  all  time  steps.  In  other  words,  for  a  deformable  mesh, 
it  is  not  suitable  to  use  the  adjoint  method  to  solve  the  inverse  problem.  The  second  problem 
is  that  the  observed  information  related  with  interaction  discharge,  such  as  the  surface  water 
stage  in  the  wetland  in  the  SV5  case,  is  difficult  to  be  used. 
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As  shown  in  Figure  4-4,  to  finish  each  iteration  of  searching  parameters,  four  steps 
including  solving  the  forward  problem  and  the  adjoint  state  equations  must  be  executed.  For 
one  parameter  in  the  case  here,  the  above  method  seems  excessive.  The  "trial  and  error" 
method  is  clearer  and  faster. 
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Figure  4-4.  Vertical  cross-section  of  the  adjoint  state  distribution  at  t=10  hours 


CHAPTER  5 
THREE-DIMENSIONAL  VARIABLY-SATURATED  GROUNDWATER  MODEL 


5.1  Introduction 

The  limitations  associated  with  saturated  subsurface  flow  models  have  inspired 
hydrologists  to  develop  variably  saturated  groundwater  models  that  consider  both  saturated 
and  unsaturated  flow  (e.g.  Voss  1984,  van  Genuchten  1980,  Kirkland  et  al.  1992,  Cox  et  al. 
1994,  and  Boufadel  et  al.  1996).  Flow  in  variably-saturated  media  has  been  investigated  by 
a  number  of  investigators  over  the  past  several  decades.  Rubin  (1968)  used  the  finite- 
difference  and  alternating-direction-implicit  (ADI)  methods  to  solve  the  two-dimensional 
Richards'  equation.  Freeze  (1971)  used  the  finite-difference  and  line  successive  over- 
relaxation  methods  to  simulate  three-dimensional  variably  saturated  flow.  Neuman  (1973) 
initiated  the  use  of  the  finite-element  method  to  model  two-dimensional  variably  saturated 
flow  in  1973.  Huyakorn  el  al.(1984,  1986)  developed  the  SATURN  two-dimensional  and 
three-dimensional  variably-saturated  flow  models  using  finite  elements.  Yeh  (1 992)  used  the 
finite-element  method  and  Gauss  elimination  (or  point  iteration)  methods  to  solve  the  three- 
dimensional  Richards'  equation.  His  model  is  the  well  known  three-dimensional 
FEMWATER  model.  Other  significant  contributions  to  the  area  of  variably-saturated  flow 
can  be  found  in  Brutsaert  (1971),  Narasimhan  (1978),  Cooley  (1983),  Kirkland  et  al.(1992), 
Clement  et  al.(1994)  and  Bates  et  al.  (2000),  Zhang  (2000),  Thorenz,  2002,  and  Russo 
(2003). 
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The  general  modeling  approach  has  consisted  of  extending  Darcy's  law  to  unsaturated 
media  and  combining  it  with  the  continuity  equation  to  obtain  Richard' s  equation.  Empirical 
relations  are  then  used  to  relate  unsaturated  hydraulic  conductivity  to  the  soil  moisture 
content  and  pressure  head  [Brooks  and  Corey  1966;  and  Van  Genuchten  1980].  Although 
there  are  many  papers  on  variably-saturated  flow  modeling,  few  have  investigated  variably- 
saturated  flow  using  deformable  finite  elements. 

This  chapter  describes  an  efficient  three-dimensional  finite-element  variably- 
saturated  groundwater  model  that  is  suitable  for  characterizing  wetland/aquifer  interactions 
using  a  deformable  mesh  and  real-time  boundary  conditions.  The  model  was  tested  using 
experimental  data  from  Vauclin  (1979)  and  then  applied  to  simulate  surface  stage  and 
groundwater  level  variations  measured  in  an  isolated  wetland  in  South  Florida.  A  comparison 
was  made  between  the  simulation  results  predicted  by  the  variably  saturated  flow  model  and 
that  of  the  saturated  flow  model  (described  in  Chapter  3)  for  the  isolated  wetland  field  study. 

5.2  Three-Dimensional  Finite-Element  Variably-Saturated  Groundwater  Model 
5.2.1  Governing  Equation 

Eq.  5-1  describes  variably-saturated  groundwater  flow,  assuming  that  the  fluid  is 
incompressible,  the  system  is  isothermal,  and  the  air  phase  is  infinitely  mobile: 


a*. 


""   'X±e=|(9A).  WW  ^ 


K«Kr»— 

)       ) 


where  S  w  is  the  degree  of  saturation,  equal  to  0/0s;  0  is  the  water  content;  0S  is  the  saturated 
water  content;  Ky  is  the  saturated  hydraulic  conductivity  tensor  [L/T];  Krw=Krw(Sw  ) 
represents  the  relative  permeability  of  the  medium,  which  is  a  function  of  the  degree  of  water 
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saturation  Sw;  t|/  is  pressure  head  [L];  and  z  is  elevation  head  [L].  In  order  to  solve  this 
nonlinear  flow  equation,  constitutive  relations  must  be  established  that  relate  the  primary 
unknowns  ijj  and  Sw.  Such  relations  are  often  given  by  a  prescribed  functional  relationship 
that  is  determined  experimentally.  The  following  functional  relation  given  by  Van  Genuchten 
(1980)  was  used  in  this  development: 


se 


1 


U(aPcy 


m  1 

,  m  =  \-—  ;     0<m<\  (5.2) 

n 


where:  a  and  m  are  parameters  obtained  from  a  fit  of  above  equations  to  experimental 
results;  Pc  is  the  capillary  pressure  ( Pc  =  -  t|j);  and: 


S  -S         0-0 

o  _     w       wr  _  r 

'    \-S         0  -0  (5'3> 

■wr  s        r 


where  S  wr  is  the  residual  saturation;  Se  is  effective  saturation;  and  0r  is  residual  moisture 
content. 

The  relative  permeability  is  obtained  from  the  effective  saturation  using  the  following 
relationship  from  Van  Genuchten  (1980): 


k    -Sm 


l-(l-5,1/m> 


(5-4) 


In  using  Eq.  5-4,  hysteresis  was  not  considered. 
5.2.2  Linearizing  the  Governing  Equations 

To  solve  the  governing  equations,  the  nonlinear  Eq.  5-1  was  linearized.  Substituting 
Eq.  5-3  into  the  right-hand-side  of  Eq.  5-1  yields: 


90 


<m)  _  ae  _  d[(Qs-er)se+er] 

dt  dt  dt 


(5-5) 


Combining  Eq.  5-5  and  Eq.  5-2  produces: 


dt 


=  (6-0> 


1 

(_1)       a»nPnldP<<SSd* 

l  +  anP" 

(5-6) 


where: 


dPc__dty  __d(ty+z) 
dt         dt  dt 


(5-7) 


Thus,  Eq.  5-5,  Eq.  5-6  and  Eq.  5-7  yield: 


w 


/,(e,,er,i»,a,pc,w,ve)^|^ 


(5-8) 


By  substituting  Eq.  5-6  into  the  left  hand  side  of  Eq.  5-1,  the  following  linearized  equation 
was  obtained: 


d  ( K  K    d(ty+z)^ 


dx. 


lj      rw 


dx. 


±Q  = 


j    > 


(9,-6> 


(  ,  \  m-\ 


\  +  a"P 


n  n  " 
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(\  +  a"Pcn) 


annP"~l+SS 

c  e    s 


d(ty+z) 


dt 


J  J  =  1,2,3 


(5-9) 
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The  governing  Eq.  5-9  was  solved  using  the  finite-element  method. 
5.2.3  Boundary  Conditions  for  Variably-Saturated  Model 

Figure  5-1  shows  a  hypothetical  wetland  consisting  of  a  surface-water  body  lined 
with  a  peat  layer,  and  beneath  the  wetland  is  an  aquifer.  The  lower  boundary  of  the  aquifer 
is  assumed  to  be  an  impermeable  unit.  As  to  the  top  boundary  condition  (IA,  BJ),  constant 
moisture  or  a  moisture  flux  can  be  specified.  In  the  case  of  rainfall,  the  net  rainfall  would 
be  a  specified  flux  as  long  as  Pprecipitation  <  Ks.  For  the  case  of  no  precipitation,  the  moisture 
flux  is  equal  to  actual  evapotranspiration. 

As  to  the  coupling  boundary  condition  controlled  by  wetland  water  level,  the  same 
way  used  in  saturated  flow  model  was  used  in  this  model.  The  initial  water  table  of  the 
wetland  and  pumping  discharge  were  taken  as  the  only  inputs.  During  pumping,  the  pumping 
discharge  lowered  the  wetland  water  level,  and  it  also  triggered  the  discharge  from  aquifer 
supplying  to  wetland  through  the  peat  layer.  In  each  simulation  time  step,  the  pumping 
discharge,  the  estimated  discharge  from  aquifer  to  wetland,  the  initial  water  table  of  wetland, 
and  the  relation  between  the  volume  of  wetland  and  wetland  water  level  were  put  together, 
and  the  water  table  at  the  end  of  the  time  step  was  estimated.  The  estimated  water  table  at 
the  end  of  the  time  step  was  again  put  back  to  the  model  of  the  variably-saturated 
groundwater  in  the  aquifer  and  peat  layer;  it  corrected  the  discharge  from  aquifer  to  wetland 
estimated  before.  Therefore,  several  iterations  in  each  time  step  were  necessary  to  get  the 
accurate  estimated  discharge  from  aquifer  to  wetland  and  the  convergent  solution  of  water 
table  in  wetland  at  the  end  of  the  time  step.  While  the  pumping  was  stopped,  the  pumping 
discharge  was  taken  as  zero.  Thus,  the  boundary  condition  at  any  time  unknown  in  the 
beginning  depended  on  the  modeling  results  of  the  coupling  problem  itself  before  this  time. 
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AB,  EG,  FH  or  CD  are  head  boundaries  T2 
GH  or  EC,  DF  are  no  flow  boundaries  T  , 
EA,  BF  are  head,  flow  boundaries  T  , 

Figure  5-1.  Sketch  showing  setup  of  the  variably  saturated 
groundwater  flow  problem  around  a  wetland 

The  determination  of  evapotranspiration  (ET)  is  a  two-step  process  in  the  model 
based  on  Fraisse  and  Campbell  (1997).  First,  the  daily  potential  evapotranspiration  (PET)  is 
calculated  in  terms  of  atmospheric  data.  Maximum  evaporation  depends  on  climatological 
factors  which  include  net  radiation,  temperature,  humidity,  and  wind  velocity.  The  PET 
represents  the  maximum  amount  of  water  that  will  leave  the  soil  system  by  evaporation  when 
there  is  a  sufficient  supply  of  soil  water.  After  PET  is  calculated,  checks  are  made  to 
determine  if  ET  is  limited  by  soil  water  conditions.  If  soil  water  conditions  are  not  limiting, 
ET  is  set  equal  to  PET.  When  PET  is  higher  than  the  amount  of  water  that  can  be  supplied 
from  the  soil  system,  ET  is  set  equal  to  the  smaller  amount.  Maximum  evaporation  is  equal 
to  potential  evaporation.  The  following  relationship  between  ET  and  PET  developed  by 
Norero  (1969)  was  used  in  the  groundwater  flow  model: 
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ET=T 


1  + 


PET 
(    i  \k 


(5-10) 


and 


2.56 


t  Wmin 


log 


* 


(5-11) 


\     ■  Wmax ) 


where  k  is  a  constant,  i|j  is  the  soil  water  potential  or  metric  pressure  in  the  root  zone,  i|/*  is 
the  value  of  i|/  when  ET=  0.5PET,  i|/Wmm  is  the  value  of  i|i  that  corresponds  to  PE/PET=0.05, 
and  tj;Wnm  is  the  value  of  i|/  that  corresponds  to  PE/PET=0.95.  Norero  (1969)  used  iJ;Wnrai  =- 
1350  joules/kg  and  t|;Wmax=-200  joules/kg.  Figure  5-2  illustrates  the  relative  evapotranspiration 
as  a  function  of  soil  water  potential. 

ET/PET 


10      8       6      4 


SOIL  WATER  POTENTTIAL,  tfr,  BARS 


Figure  5-2.  Schematic  of  relative  evapotranspiration  (ET/PET), 
as  affected  by  soil  water  potential,  i|/  [Norero  1969]. 
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5.3  Finite-Element  Method  for  Modeling  Three-Dimensional  Variably-Saturated  Flow 

Given  the  complexity  of  the  unsaturated  flow  Eq.  5-9,  it  is  not  readily  seen  that  it  has 
the  same  form  as  the  governing  equation  for  satuarted  flow  (Eq.  3-1).  However,  if 
f(0s,  0r,  m,  a,  pc,  n)  from  Eq.  5-8  is  used  in  lieu  of  specific  storativity  Ss  in  the  saturated 
groundwater  equation,  the  form  of  both  the  saturated  and  unsaturated  flow  equations 
becomes  the  same.  As  a  result,  the  strategy  taken  to  develop  an  finite-element  model  for 
saturated  flow  (Chapter  3)  was  applied  again  to  develop  an  appropriate  equation  for  an  finite- 
element  model  for  variably  saturated  flow.  Similar  to  what  was  done  in  Chapter  3,  assuming 
that  the  basis  function  w  is  chosen  from  V  space,  the  follow  equation  is  always  true: 


w[v-(KKrwVh)]dQ±  I  wQdQ 


=  \  wft(Qs,Qr,m,a,Pc,n)-±dQ,    zj-1,2,3  (5.12) 


* 


Manipulation  of  the  first  term  in  the  left  hand  side  of  Eq.  5-12  begins  with  the  application 
of  Green's  theorem,  which  is  an  expanded  form  of  the  divergence  theorem.  This  converts  the 
integral  into  two  terms,  one  of  which  is  evaluated  only  at  the  surface  of  the  region  to  be 
simulated.  Application  of  the  Galerkin  criterion  and  Green's  theorem  to  Eq.  5-12  leads  to 
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w[(AXwV/*)-fl]rfr+|    w[(£X^VA)-ii]rfr-l  {'7w)-(KKnvVh)dQ 
r,  Jr2  Jo. 


I    w0rfQ  =  |   W^(0i50r,7H,a,i>)^</Q  (5-13) 


The  first  term,  which  depends  on  the  moisture  flux  on  the  surface,  is  no  longer  zero  as  it  is 
in  the  saturated  model.  For  head  boundaries  that  are  used  along  T2,  w  is  equal  to  zero, 
because  it  belongs  to  V  space,  Thus,  with  the  second  term  set  equal  to  zero,  Eq.  5-13  can  be 
rewritten  as: 

w[(JTJ^VA)-||]rfr-J  {Vw)-(KKnvVh)dQ±l  wQdQ 


wft(e,e„m,a,P(i,ri&dQ  (5-14) 


An  approximate  solution  for  h(x,  y,  z,  t)  can  be  assumed  to  have  the  following  form: 


WW 


=X>/'W,*K    ,  NteV  (5-15) 


and: 
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1  =  1 


(5-16) 


where  h(e)  is  the  approximate  solution  for  hydraulic  head  within  element  e,  N/e)  are 
interpolation  functions  for  each  node  within  element  e,  p  is  the  number  of  nodes  within 
element  e,  and  d;(t)  are  the  unknown  values  of  hydraulic  head  for  each  node  within  element 
e  and  at  time  t.  At  the  head  boundaries,  I\,  d;(t)  =  h(x,y,z)|r2.  w/e)  is  the  weighting  function 
for  node  i  and  the  limits  of  integration  are  chosen  to  represent  the  volume  of  the  elements. 
In  Galerkin's  method,  the  weighting  function  for  each  node  in  the  element  is  set  equal  to  the 
interpolation  function  for  the  node,  Wj(e,=  N/e). 

When  the  approximate  solutions  for  Eq.  5-15  and  Eq.  5-16  are  substituted  into  Eq. 
5-14,  Eq.  5-17  is  obtained  which  is  not  satisfied  exactly,  and  an  error  or  residual  occurs  at 
every  point  in  the  problem  domain.  If  it  is  also  assumed  that  the  aquifer  is  isotropic,  that  is 
to  say,  values  of  saturated  hydraulic  conductivity  in  the  three  coordinate  directions  are 
constant  within  an  element  (but  can  vary  from  one  element  to  the  next),  the  contribution  of 
any  element  e  to  the  residual  at  a  node  i  associated  with  the  element  can  be  described  through 
the  flow  system  of  Eq.  5-17: 


R?\t) 


*f(0 


K^K 


4e\t) 


4%) 


(5-17) 


where: 
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Combining  Eq.  5-17-Eq.  5-20  from  each  element  together,  global  equations  are  obtained: 


98 


\KK 


dx(t) 

♦M 

dt 

-M- 

^(0 

dp(t) 

ddp«) 

dt 

*,(') 

(5-21) 


By  setting  the  residuals  equal  to  zero,  the  final  system  of  ordinary  differential  equations  is 
obtained: 


K] 


dx(t) 

3/ 

0 

+  [c] 

-M-o 

dp(t) 

0 

(5-22) 


The  system  of  Eq.  5-22  representing  the  variably  saturated  flow  system  was  solved  using 
finite  differences. 

5.4  Transient,  Variably-Saturated  Water-Table  Recharge  Example  Test 

An  example  problem  was  formulated  using  laboratory  data  to  verify  the  performance 
of  the  model.  Details  concerning  the  following  laboratory  data  and  the  experiment  were 
presented  by  Vauclin  et  al.  (1979).  The  flow  domain  consisted  of  a  rectangular  soil  slab, 
6.00  m  by  2.00  m  (Figure  5-3),  with  an  initial  horizontal  water  table  located  at  a  height  of 
0.65  m.  At  the  soil  surface,  a  constant  flux  of  q  =  3.55  m/day  was  applied  over  an  infinite 
zone  1  m  wide.  The  remaining  soil  surface  was  covered  to  prevent  evaporation  losses. 
Because  of  symmetry,  it  was  necessary  to  model  only  one  side  (for  this  test  the  right  side)  of 
the  flow  domain.  No-flow  boundaries  were  set  on  the  bottom  and  along  the  axis  of  symmetry 
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(  Figure  5-3).  At  the  soil  surface,  the  constant  flux  of  3.55  m/day  was  applied  over  the  left 
0.50  m  of  the  top  of  the  modeled  domain.  The  remaining  soil  surface  represented  a  no-flow 
boundary  condition.  A  no-flow  boundary  was  used  above  the  water  table  on  the  right  side  of 
the  domain  (Figure  5-3). 


Infiltration  Zone 


^■jr 1  j  ^/So\X  Sufac 
•A    B  C~ 


Initial  Position  of  the  Water  Table 

J?  E 


Lateral 


I  Trench 


///////////////////////////////// 
Impervious  Substratum 

DE  is  head  boundary 

AB  is  flux  boundary 
BC,CD,AF,FE  are  no  flux  boundary 

Figure  5-3.  Schematic  diagram  of  the  flow  domain 
From  Vauclin  et  al.  (1979),  the  characteristics  of  the  experimental  soil  properties 
included  a  saturated  hydraulic  conductivity  of  8.40  m/day,  a  porosity  of  r\  =  0S  =  0.30,  and 
a  residual  saturation  of  0  =  0.01.  When  Eq.  5-2  was  fitted  to  the  water  retention  and  the 
relative  hydraulic  conductivity  data  given  by  Vauclin  et  al.  (1979),  estimated  values  for 
model  parameters  a  and  n  were  a  =  3.3  m'1  and  n  =  4.1.  Specific  storage  was  neglected  for 
this  problem  because  changes  in  storage  are  facilitated  by  the  filling  of  pores,  which 
overshadows  the  effects  of  compressibility  for  the  vertical  extent  of  this  flow  domain.  Hence, 
the  specific  storage  coefficient  was  set  to  zero. 

The  mesh  spacing  used  in  the  test  problem  in  the  horizontal  x-  and  vertical 
z-directions  was  Ax  =  0.10  m  and  Az  =  0.05m,  respectively.  The  soil  system  was  assumed 
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to  be  initially  at  hydrostatic  equilibrium  with  respect  to  the  water  table  throughout  the  flow 
domain.  Transient  simulations  were  performed  using  time  steps  that  varied  from  At  =  180 
secondes  to  At  =  360  secondes. 

Transient  positions  of  the  water  table  are  compared  with  the  experimental  results 
presented  by  Vauclin  et  al.  (1979)  in  Figure  5-4  and  Figure  5-5.  These  figures  show  excellent 
agreement  between  the  transient  water-table  positions  predicted  by  the  variably  saturated 
model  and  those  obtained  from  Vauclin  et  al.  (1979). 
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Figure  5-4.  Water  moisture  distribution  (x  :  y  =  2  : 1)  at  4  hours 
(Observed  points  are  from  Vauclin  et  al.1979) 
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Figure  5-5.  Water  moisture  distribution  at  8  hours  (x  :  y  =  2:1) 
(Observed  points  are  from  Vauclin  et  al.1979) 


5.5  Application  on  SV5  Wetland 
5.5.1  Finite-Element  Mesh  for  Three-Dimensional  Variably-Saturated  Flow  Model 

Unlike  the  saturated  groundwater  flow  model  (Chapter  3),  matric  pressure  and  water 
content  may  vary  rapidly  in  time  and  space  in  the  variably-saturated  model,  which  produces 
rapid  spatial  and  temporal  variation  in  the  variably-saturated  flow  hydraulic  conductivity. 
In  order  to  simulate  the  moisture  or  pressure  variation  accurately,  a  fine  numerical  mesh  must 
be  used  over  zones  of  rapid  change  to  determine  flow  parameter  variations  accurately.  Like 
the  saturated  flow  model,  the  SV5  (Wise  et  al.  2000)  simulation  domain  was  divided  into  two 
main  areas.  One  area  was  the  aquifer;  the  other  area  was  the  wetland  peat  layer.  Two 
parameters,  kaquifer  and  k^,  were  used  to  represent  their  hydraulic  conductivities.  The 
numerical  mesh  was  26x21x3.  The  top  10  layers  of  mesh  represent  about  top  1.55  m  of  the 
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aquifer.  Among  the  top  10  layers,  the  upper  4-5  layers  were  given  a  finer  descretization,  and 
the  remaining  layers  of  the  top  10  were  given  a  coarse  descretization.  The  bottom  10  layers 
represented  about  15  meters  of  saturated  aquifer.  The  mesh  in  the  saturated  portion  of  the 
flow  field  was  considerably  coarser. 


-    15 


-    10 


N 


Figure  5-6.  Finite-element  mesh  of  the  wedge  part  for  variably  saturated  flow  model 

(Unit  is  m) 

Unlike  the  model  for  saturated  flow,  certain  of  the  boundaries  surrounding  the 

simulation  area  are  fixed.  The  top  boundaries  are  the  aquifer  surface  or  the  wetland  bottom. 

The  mesh  deforms  during  a  simulation  because  the  elevation  of  the  surface  water  changes. 

The  model  seeks  to  locate  the  same  elevation  as  the  surface-water  stage  and  at  a  location 

where  a  horizontal  line  defines  where  the  surface-water  surface  intersects  with  the  bottom 

of  the  wetland.  As  the  mesh  deforms  during  a  simulation,  a  subroutine  program  is  used  to 

determine  whether  a  cell  is  in  the  peat  layer  or  the  underlying  aquifer  and  then  the 

corresponding  hydraulic  conductivity  and  other  parameters  of  the  element  are  assigned.  As 

the  mesh  deforms,  some  elements  are  enlarged  and  some  of  them  are  compressed,  but  the 
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total  area  of  the  model  domain  does  not  change. 
5.5.2  Modeling  Results  and  Summary 

For  the  SV5  wetland  simulation,  the  same  hydraulic  conductivity  values  used  in  the 
saturated  model  were  used  in  the  variably  saturated  flow  model.  Thus,  Kxx=Kyy=8.5  m/day; 
K„=8.5  m/day;  the  anisotropy  ratio  Kh/Kv=  1.0  in  aquifer;  and  Kp=  0.2  m/day  in  the  peat 
layer  were  used  throughout  the  model  domain.  The  specific  storativity  of  Ss  =0.00001  (1/m) 
and  the  saturated  water  content  and  residual  water  content  in  the  peat  layer  0S  =  0.7;  0r  =  0.6; 
n=5.0;  m=0.8;  a  =  0.001  (1/cm)  =  0.1  (1/m);  and  in  the  aquifer  0S  =  0.34;  0r  =  0.18;  n=9.1; 
m=0.89;  a  =  0.01(l/cm)  =  1  (1/m)  were  also  used,  based  on  Walser  et  al.  (1998)  and  Wise 
et  al.  (2000).  Figure  5-7  shows  the  simulated  piezometric  head  variation  at  selected  depths 
beneath  the  wetland. 

When  surface  water  was  pumped  from  the  wetland,  the  water  stage  decreased,  and 
the  associated  piezometric  head  changes  were  greatest  immediately  beneath  the  wetland.  The 
time  needed  for  the  wetland  to  recover  from  the  stress  was  much  longer.  Figure  5-8  shows 
that  the  aquifer  discharge  to  the  wetland  increased  as  the  surface-water  component  of  the 
wetland  was  drained  and  that  the  aquifer  discharge  decreased  a  little  after  the  pumping 
stopped.  Figure  5-9  shows  the  transient  variation  of  the  inundated  wetland  area.  During 
pumping,  the  inundated  area  decreased  from  2,800  m2  to  700  m2.  During  the  recovery  period, 
the  inundated  area  increased  slowly.  Figure  5-10  shows  that  the  model  was  able  to  fit  the 
observed  aquifer  well,  which  is  the  SFWMD  well  screened  at  a  depth  of  six  meters  and 
grouted  from  the  surface  to  a  depth  of  five  meters.  In  a  comparison  prediction,  results  from 
both  the  saturated  model  and  the  variably  saturated  flow  model  are  consistent  (Figure  5-1 1). 
The  piezometric  head  distribution  beneath  the  wetland  is  shown  in  Figure  5-12  and  Figure 


104 


5-13.  Tables  5-1  through  5-5  show  the  simulated  heads,  the  corresponding  areas,  and  the 
distances  which  were  used  to  calculate  the  exchange  discharge. 
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Figure  5-7.  Calculated  piezometric  head  beneath  the  wetland  using 
the  three-dimensional  variably  saturated  groundwater  flow  model 
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Figure  5-8.  Model  predicted  discharge  from  the  aquifer  to  the  overlying  wetland 
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Figure  5-9.  Inundated  area  of  wetland  during  simulation 
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Figure  5-10.  Simulated  and  observed  heads  in  the  wetland  and  in  aquifer 


106 


4.25 

4.2 

£     4.15 

I        4.1 

o 

jB      4.05 

o 

N 

o)  4 

a. 

3.95 

3.9 


\ 

S-M  (WETL) 
V-M  (WETL) 
S-M  (GW) 
V-S  (GW) 


20 


40  60 

Time  (hour) 


80 


Figure  5-11.  Simulated  heads  from  the  saturated  flow  model 
and  the  variably  saturated  flow  model 
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Figure  5-12.  Horizontal  view  of  the  three-dimensional  finite-element  mesh  and 
piezometric  head  distribution  from  the  variably  saturated  flow  model  at  a  time  of  9  hours. 
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Figure  5-13.  Enlarged  top  and  middle  part  of  the  mesh 
and  head  contour  at  a  time  of  9  hours. 


Table  5-1.  Data  for  calculating  the  discharge  to  the  wetland  from  the  aquifer  at  10  hours 


Elements 
of  Peat  Layer 

Area 

(m2) 

Head  in  Peat 
Layer   (m) 

Head  on  Peat 
Surface    (m) 

Distance  of  Two 
Nodes       (m) 

1 

6.82502 

15.7300 

15.7203 

0.04715 

2 

20.4751 

15.7301 

15.7203 

0.04718 

3 

34.1251 

15.7301 

15.7203 

0.04721 

4 

47.7747 

15.7291 

15.7203 

0.04723 

5 

61.4246 

15.7285 

15.7203 

0.04726 

6 

75.0755 

15.7289 

15.7203 

0.04729 

7 

88.7265 

15.7292 

15.7203 

0.04732 

8 

102.378 

15.7296 

15.7203 

0.04736 

9 

116.030 

15.7301 

15.7203 

0.04741 

10 

129.683 

15.7308 

15.7203 

0.04747 

Qto  wetland  at  10  hr 

1.6954  mVhr 
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Table  5-2.  Data  for  calculating  the  discharge  to  the  wetland  from  the  aquifer  at  20  hours 


Elements 
of  Peat  Layer 

Area 
(m2) 

Head  in  Peat 
Layer   (m) 

Head  on  Peat 
Surface    (m) 

Distance  of  Two 
Nodes      (m) 

1 

8.39429 

15.7513 

15.7431 

0.04715 

2 

25.1829 

15.7514 

15.7431 

0.04718 

3 

41.9713 

15.7515 

15.7431 

0.04721 

4 

58.7594 

15.7506 

15.7431 

0.04724 

5 

75.5484 

15.7501 

15.7431 

0.04727 

6 

92.3379 

15.7505 

15.7431 

0.04731 

7 

109.128 

15.7508 

15.7431 

0.04735 

8 

125.919 

15.7512 

15.7431 

0.04740 

9 

142.711 

15.7517 

15.7431 

0.04746 

10 

159.504 

15.7524 

15.7431 

0.04753 

Q  to  wetland  at  20hr 

1.8041  mVhr 

Table  5-3.  Data  for  calculating  the  discharge  to  the  wetland  from  the  aquifer  at  40  hours 


Elements 
of  Peat  Layer 

Area 
(m2) 

Head  in  Peat 

Layer   (m) 

Head  on  Peat 
Surface    (m) 

Distance  of  Two 
Nodes       (m) 

1 

11.0360 

15.7876 

15.7813 

0.04715 

2 

33.1081 

15.7876 

15.7813 

0.04719 

3 

55.1798 

15.7870 

15.7813 

0.04722 

4 

77.2516 

15.7866 

15.7813 

0.04726 

5 

99.3248 

15.7869 

15.7813 

0.04729 

6 

121.399 

15.7871 

15.7813 

0.04734 

7 

143.474 

15.7874 

15.7813 

0.04740 

8 

165.550 

15.7877 

15.7813 

0.04747 

9 

187.628 

15.7882 

15.7813 

0.04754 

10 

209.707 

15.7889 

15.7813 

0.04763 

Q  to  wetland  at  40hr 

1.8693  mVhr 
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Table  5-4.  Data  for  calculating  the  discharge  to  the  wetland  from  the  aquifer  at  60  hours 


Elements 
of  Peat  Layer 

Area 
(m2) 

Head  in  Peat 
Layer  (m) 

Head  on  Peat 
Surface    (m) 

Distance  of  Two 
Nodes      (m) 

1 

13.1576 

15.8171 

15.8120 

0.04715 

2 

39.4727 

15.8171 

15.8120 

0.04719 

3 

65.7872 

15.8166 

15.8120 

0.04723 

4 

92.1052 

15.8163 

15.8120 

0.04723 

5 

118.419 

15.8165 

15.8120 

0.04727 

6 

144.737 

15.8167 

15.8120 

0.04731 

7 

171.057 

15.8170 

15.8120 

0.04737 

8 

197.378 

15.8173 

15.8120 

0.04744 

9 

223.701 

15.8177 

15.8120 

0.04752 

10 

250.027 

15.8184 

15.8120 

0.04772 

Q  to  wetland  at  60hr 

1.8396  m7hr 

Table  5-5.  Data  for  calculating  the  discharge  to  the  wetland  from  the  aquifer  at  80  hours 


Elements 
of  Peat  Layer 

Area 
(m2) 

Head  in  Peat 
Layer    (m) 

Head  on  Peat 
Surface    (m) 

Distance  of  Two 
Nodes       (m) 

1 

14.9321 

15.8419 

15.8377 

0.04715 

2 

44.7963 

15.8419 

15.8377 

0.04720 

3 

74.6597 

15.8415 

15.8377 

0.04724 

4 

104.525 

15.8412 

15.8377 

0.04728 

5 

134.391 

15.8415 

15.8377 

0.04733 

6 

164.259 

15.8417 

15.8377 

0.04739 

7 

194.129 

15.8419 

15.8377 

0.04747 

8 

224.001 

15.8422 

15.8377 

0.04756 

9 

253.875 

15.8426 

15.8377 

0.04767 

10 

283.753 

15.8432 

15.8377 

0.04779 

Q  to  wetland  at  80hr 

1.7653  mVhr 

110 

The  saturated  model  and  the  variably  saturated  model  are  different  in  several  aspects. 
First  of  all,  the  variably  saturated  flow  model  requires  that  the  moisture  flux  be  specified  on 
the  soil  surface  as  the  boundary  condition.  This  moisture  flux  on  the  soil  surface  induces  the 
head  distribution  in  the  vadose  zone  that  is  not  described  in  the  saturated  flow  model. 
Second,  the  empirical  formula  used  to  characterize  the  matric  pressure  as  a  function  of  the 
degree  of  saturation  may  not  be  accurate,  and,  as  a  result,  it  displaces  the  phreatic  surface 
(where  pressure  head  is  zero)  to  a  location  not  consistent  with  observations. 

Table  5-6  shows  the  sensitivities  of  the  heads  related  to  charges  in  its  hydraulic 
conductivities  in  the  peat  layer  and  aquifer  and  to  changes  in  the  specific  storativity.  The 
hydraulic  conductivities  in  both  the  peat  layer  and  in  the  aquifer  are  important  in  the  change 
of  head  distribution.  Changes  in  the  specific  storativity  do  not  affect  the  head  very  much 
relatively . 

Table  5-6  Sensitivity  calculated  by  variably  saturated  groundwater  flow  model 


Parameter  (PR) 

Percentage  varied 

APR 

II  h,-h2ll  (m) 

llh,-h2ll/APR 

Ka(m/day) 

-0.1 

-0.85 

0.5700E-02 

-0.6706E-02 

Ka(m/day) 

+0.1 

+0.85 

0.5090E-02 

0.5988E-02 

Kptm/day) 

-0.1 

-0.02 

0.5155E-02 

-0.2578 

KpCm/day) 

+0.1 

+0.02 

0.4665E-02 

0.2333 

S0 

-0.5 

-0.000005 

0.2990E-04 

-5.980 

So 

+0.5 

+0.000005 

0.1998E-04 

3.996 

CHAPTER  6 
SUMMARY  AND  CONCLUSIONS 


This  dissertation  describes  an  analytical  formula  for  groundwater  variation  affected 
by  a  flood  wave  in  a  wetland  and  also  describes  the  development  of  two  three-dimensional 
finite-element  method  numerical  models,  one  of  which  is  a  saturated  groundwater  flow 
model  and  the  other  a  variably  saturated  groundwater  flow  model.  Both  models  use  a 
deformable  mesh  and  real-time  coupling  for  the  wetland-aquifer  boundary  along  the  interface 
between  the  surface  water  and  groundwater.  In  the  saturated  groundwater  flow  model,  a 
predictor-corrector  method  is  used  to  find  the  real-time  boundary  condition  that  satisfies  both 
the  head  condition  and  flux  condition  at  the  same  time.  The  saturated  groundwater  flow 
model  was  verified  using  two  analytical  solutions;  one  is  the  Theis  solution,  and  the  other 
is  the  analytical  solution  found  in  Chapter  2  of  the  dissertation,  which  describes  groundwater 
flow  around  a  circular  wetland.  The  saturated  groundwater  flow  model  was  further  used  to 
simulate  the  interaction  problem  of  groundwater  and  the  SV5  wetland  in  south  Florida. 

The  three-dimensional  finite-element  variably-saturated  flow  model  was  developed 
based  on  the  linearized  Richard  equation.  This  model  is  different  from  previous  variably- 
saturated  flow  models  in  the  deformable  three-dimension  finite-element  mesh  and 
dynamically  interacting  boundary  conditions  applied  in  this  model.  The  predictor-corrector 
method  is  also  used  in  the  variably-saturated  groundwater  flow  model  to  find  the  real-time 
boundary  conditions  that  satisfies  both  the  head  condition  and  flux  condition  at  the  same 
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time.  The  variably-saturated  groundwater  flow  model  was  verified  based  on  the  field  and 
laboratory  work  described  by  Vauclin  et  al.  (1979). 

The  following  specific  conclusions  can  be  drawn  from  this  dissertation: 

1.  The  analytical  model  describing  the  interaction  between  a  flood  wave  in  a  lake  and 

surrounding  groundwater  is  an  original  development  that  can  be  used  to  simulate  the 
groundwater  fluctuations  caused  by  a  transient  wetland  surface  water  stage  regardless 
of  the  function  used  to  characterize  the  transient  variation  of  the  wetland  stage. 

2.  Considering  that  the  area  of  the  water  surface  changes  quickly  because  most  wetlands 

have  a  shallow  bank  slope,  the  saturated  groundwater  flow  model  has  the  capability 
of  successfully  tracking  the  interaction  boundary  between  groundwater  and  a 
wetland.  The  pumping  from  the  SV5  wetland  and  the  recovery  of  the  wetland  water 
surface  afterwards  were  closely  simulated  by  this  model,  and  thus  the  interacting 
boundaries  were  effectively  traced  by  the  deformable  mesh. 

3.  The  variably  saturated  groundwater  flow  model  also  successfully  simulated  the 

groundwater  flow  affected  by  the  pumping  from  the  S  V5  wetland  and  the  water  level 
recovery  in  the  wetland  after  pumping.  The  rapidly  moving  boundary,  i.e.,  the  special 
phenomenon  that  occurs  in  the  groundwater  and  wetland  interaction  problem,  is 
effectively  dealt  with  by  this  model  also  using  a  deformable  mesh. 
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